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1.0  INTRODUCTION 


Diffusers  are  an  important  component  in  many  areas  of  application, 
such  as  propulsion  systems,  wind  tunnels,  test  facilities,  etc.  Yet  the 
prediction  of  diffuser  flows  remains  one  of  the  most  difficult  fluid  dy- 
namics problems,  especially  when  the  inlet  conditions  to  the  diffuser 
are  highly  nonuniform.  Since  the  optimum  operating  condition  for  a 
diffuser,  i.e.  , maximum  pressure  recovery,  has  been  shown  experi- 
mentally to  occur  with  some  flow  separation  (Ref.  1),  a realistic  solu- 
tion for  the  flow  field  can  be  obtained  only  by  solving  the  full  Navier- 
Stokes  equations.  In  addition,  the  diffuser  flows  of  practical  interest 
are  turbulent  in  nature.  The  performance  of  a diffuser  depends  not 
only  on  the  shape  of  the  inlet  velocity  profile  but  also  on  the  turbulence 
level.  The  modeling  of  turbulence  with  the  large  pressure  gradient 
existing  in  the  diffuser  flows  requires  a more  sophisticated  approach 
than  the  use  of  simple  eddy  viscosity  models  (Ref.  2). 

Currently,  the  diffuser  design  information  is  obtained  almost  solely 
from  empirical  data  (Refs.  1 and  3).  Many  of  the  available  diffuser 
performance  maps  and  correlations  provide  only  static  pressure  re- 
covery. Very  few  detailed  turbulence  properties  of  diffuser  flows,  es- 
tablished experimentally,  are  available  (Refs.  2 and  4)  for  non-separated 
cases.  Data  for  separated  diffuser  flows  is  even  more  sparse  because 
conventional  instruments,  such  as  the  hot  wire  anemometer,  cannot  pro- 
vide meaningful  measurement  in  regions  where  the  flow  direction  re- 
verses with  time.  This  situation  may  improve  in  the  future  as  the  re- 
cently developed  laser  velocimeter  (LV)  becomes  more  available  and 
reliable  (Refs.  5 and  6)  so  that  the  flow  field  data  can  be  obtained  for 
use  in  the  development  and  verification  of  analytical  prediction  methods. 

In  the  past,  the  flow  in  a diffuser  has  been  analyzed  by  assuming 
that  diffuser  flow  can  be  approximated  by  a thin  boundary  layer  adjacent 
to  the  wall  and  an  inviscid  core  in  the  center  of  the  diffuser.  The 
boundary-layer  equation  and  the  inviscid  core  equation  are  then  solved 
with  or  without  interaction  (Refs.  7,  8,  and  9).  No  rigorous  method  is 
available  to  analyze  diffuser  flows  with  a highly  nonuniform  inlet  pro- 
file with  or  without  separation  (Refs.  10  and  11). 

The  purpose  of  the  investigation  reported  herein  is  to  develop  nu- 
merical prediction  methods  for  the  calculation  of  turbulent,  incompres- 
sible, separated,  subsonic  diffuser  flows  with  nonuniform  inlet  conditions. 
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The  theory  development,  the  turbulence  models,  the  coordinate  trans- 
formation, and  the  numerical  finite  difference  solution  procedures  are 
presented  along  with  comparison  of  the  results  with  available  experi- 
mental data. 


2.0  GOVERNING  EQUATIONS 


2.1  NAVIER-STOKES  EQUATIONS  AND  REYNOLDS  STRESSES 


The  basic  equations  which  describe  the  motion  of  laminar  or  tur- 
bulent flow  of  an  incompressible  fluid  are  the  Navier-Stokes  equations. 
In  turbulent  flow,  it  becomes  necessary  to  use  some  averaging  proce- 
dure (or  statistical  method)  to  provide  useful  information  about  the 
gross  features  of  the  flows  which  are  random  in  nature.  Among  the 
•methods  available,  the  time -averaged  method  has  been  widely  used  for 
constant  density  non -reacting  flows.  The  resultant  equations,  as  shown 
by  Osborne  Reynolds  (Ref.  12),  can  be  written  in  2-D  Cartesian  or 
cylindrical  coordinates: 


Continuity  Equation 


-&-('**) + fjr<rSv) 


(1) 


X- Momentum  Equation 


* lf)l} **>} 


(2) 


Y -Momentum  Equation 


£[*  <fMf )]  ♦ srW-  * f] 


(3) 


where  if  6 = 0,  r represents  y in  Cartesian  coordinates,  and  if 
6 = 1,  r represents  the  radial  coordinate. 
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In  Eqs.  (1)  through  (3),  the  flow  variables  (u,  v,  p,  and  p)  are  time- 
averaged  quantities.  The  turbulence  quantities,  such  as  u'2,  v'2, 
w'2,  and  u'v'  are  usually  called  Reynolds  stresses. 


2.2  THE  EDDY  VISCOSITY  CONCEPT 


The  eddy  viscosity  concept,  which  relates  the  Reynolds  stresses  to 
the  product  of  the  time -averaged  velocity  gradient  and  an  eddy  viscosity, 
can  be  attributed  to  Boussinesq  (Ref.  12).  With  this  concept,  Reynolds 
stresses  can  be  defined  as 


u'*  = +Tk 

V*  -*H#  +tk 


w#2 

U'V' 

k 


. 2 y X.  4 JL  u 

2 yt  r 3 


ir 


( u'*  + v'2  + w'2  )/z 


(4) 


where  v ^ is  the  eddy  viscosity  and  k is  the  turbulent  kinetic  energy  (TKE). 


By  substituting  Eq.  (4)  into  Eqs.  (2)  and  (3),  one  obtained  the  mo- 
mentum equations  in  the  following  forms: 


U»1  + 

"ax  f»r 


(5) 


(6) 


Equations  (5)  and  (6)  are  similar  in  form  to  the  original  Navier- 
Stokes  equations  if  one  replaces  (p/p  + 2/3  k)  by  (p/p)  and  ( v + v t)  by  v. 
Thus,  methods  developed  for  solving  the  Navier-Stokes  equations  can  be 
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Used  to  solve  Eqs.  (5)  and  (6),  which  allows  both  laminar  and  turbulent 
flows  to  be  solved  within  one  numerical  framework. 


2.3  VORTICITY-STREAM  FUNCTION  FORMULATION 


Equations  (5)  and  (6)  are  coupled  nonlinear  partial  differential  equa- 
tions. It  can  be  seen  that  the  pressure  term  appears  only  in  the  gra- 
dient form,  i.  e. , 9/8x(p/p)  and  a/ar(p/p).  In  most  cases,  the  pres- 
sure distribution  is  unknown.  Therefore,  it  is  advantageous  to  remove 
the  explicit  pressure  gradient  terms  from  Eqs.  (5)  and  (6).  This  is 
done  by  a cross -differentiation  (or  a curl  operation).  The  resultant 
equation  written  in  terms  of  the  vorticity  becomes 


= 0 


<v) 


where  Q is  the  vorticity  defined  as 


Cl  a 


ra 

ir 


(8) 


The  continuity  equation  (Eq.  (1))  must  be  modified  to  obtain  a solu- 
tion because  it  is  a first-order  differential  equation  as  opposed  to  the 
vorticity  equation  (Eq.  (7)),  which  is  a second-order  equation.  A 
second-order  differential  equation  can  be  formulated  to  replace  the  con- 
tinuity equation  by  introducing  a stream  function. 


u 

V 


(9) 


The  stream  function  defined  in  Eq.  (9)  satisfied  continuity  Eq.  (1)  auto- 
matically. By  combining  Eqs.  (9)  and  (8),  a single  second-order  differ- 
ential equation  for  the  stream  function  ($ ) is  obtained 
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( 


aV 

•ax* 


' r J9r 


£ 

+ r £i=o 


(10) 


Equations  (7)  and  (10)  replace  Eqs.  (5),  (6),  and  (1)  to  form  a "Vorticity- 
Stream  Function"  formulation. 


2.4  THE  PRESSURE  EQUATION 


Since  the  pressure  term  was  eliminated  in  the  first  stage  of  the 
formulation,  it  needs  to  be  recovered  after  the  and  0 solutions  are  ob- 
tained. Two  methods  are  available,  namely:  (1)  the  integration  of  mo- 
mentum Eqs.  (5)  and  (6),  and  (2)  the  solution  of  a pressure  equation. 

The  first  method  is  straightforward.  Since  the  pressure  appeared  in 
the  momentum  equation  in  the  gradient  form,  i.  e. , did  x(p / p)  and 
3/9r(p Ip),  one  can  integrate  the  momentum  equation  from  a reference 
point  where  the  pressure  is  known  to  provide  a continuous  distribution 
along  any  path,  i.  e.  , 


X-Integration 


(11) 


r -Integration 


(12) 


The  accuracy  of  the  pressure  calculated  by  Eqs.  (11)  and  (12)  generally 
depends  on  the  numerical  integration  quadrature  in  addition  to  the  accu- 
racy of  the  numerically  calculated  variables,  such  as  u,  v,  k,  etc.  The 
method  is  commonly  used  to  evaluate  the  wall  pressure  distribution 
where  Eqs.  (11)  and  (12)  can  be  simplified  because  of  the  non-slip  wall 
boundary  condition,  i.  e.  , u = v = 0 at  the  wall.  The  approach  can  also 
provide  an  initial  pressure  distribution  to  be  used  as  an  initial  guess  in 
the  solution  of  the  pressure  equation  given  below. 
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The  second  method  requires  the  solution  of  a pressure  equation 
which  can  be  derived  by  a further  differentiation  (or  a gradient  opera- 
tion) of  momentum  Eqs.  (5)  and  (6), 


(13) 


Although  Eq.  (13)  can  provide  uniquely  defined  pressure  distributions 
as  opposed  to  the  path  dependent  integration  procedure,  one  needs  to 
solve  an  additional  equation  along  with  the  related  boundary  conditions. 
Even  with  the  pressure  distribution  uniquely  determined,  the  accuracy 
still  depends  on  the  numerical  method  used.  However,  in  contrast  to 
the  first  method,  the  error  does  not  accumulate  along  a particular  path 
of  integration. 

The  pressure  equation  (Eq.  (13))  is  a Poisson's  equation  with 
source  terms  appearing  on  the  right-hand  side  of  the  equation.  In  the 
present  analysis,  Eq.  (13)  is  solved  iteratively  by  a point  iteration 
method.  The  specification  of  the  boundary  condition  for  the  pressure 
equation  requires  special  attention  and  is  discussed  in  detail  in  Section 
6.1. 


3.0  MATHEMATICAL  MODELS  OF  TURBULENCE 


The  concept  of  the  eddy  viscosity  (i/j.)  introduced  in  the  previous 
section  simplified  the  modeling  of  the  Reynolds  stresses.  Thus,  in- 
stead of  considering  the  Reynolds  stresses  directly,  one  can  work  with 
a single  eddy  viscosity.  The  modeling  of  turbulence,  in  this  case,  is 
directly  related  to  the  modeling  of  the  eddy  viscosity.  Because  the  eddy 
viscosity  is  not  a physical  property  of  the  fluid,  one  needs  a mathemati- 
cal relation  to  calculate  it.  In  general,  the  models  of  the  eddy  viscosity 
can  be  divided  into  three  groups  depending  on  the  degree  of  sophistica- 
tion, namely,  (1)  constant  values,  (2)  algebraic  relations  to  the  mean 
flow  property  and  (3)  differential  equation  representations  of  the  turbu- 
lence properties  (Refs.  13  and  14). 
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3.1  CONSTANT  VISCOSITY  MODEL 


The  simplest  eddy  viscosity  model  is 

>>  a constant  (14) 

in  which  the  eddy  viscosity  is  given  a constant  value  throughout' the  flow 
field.  This  modeling,  although  crude,  can  provide  approximate  solu- 
tions when  used  in  the  calculation  of  certain  simple  flow  problems,  such 
as  the  far  field  of  a free  turbulent  round  jet.  The  model  is  also  useful 
in  the  early  stage  of  the  development  of  a particular  numerical  method 
or  procedure  because  the  resultant  governing  equations  are  the  same 
as  those  derived  from  laminar  flows  with  a constant  molecular  viscosity. 
The  only  difference  is  that  the  eddy  viscosity  is  usually  several  orders 
of  magnitude  larger  than  the  molecular  viscosity.  Without  the  spatial 
variation  of  the  eddy  viscosity,  one  can  concentrate  directly  on  the  sta- 
bility or  the  convergence  of  a particular  numerical  method  used  in  the 
analysis . 


The  vorticity  equation  (Eq.  (7))  in  the  present  vorticity-stream 
function  formulation  can  be  simplified  to  the  following  form  by  using 
the  constant  eddy  viscosity  model: 


(15) 


Equations  (14),  (15),  and  (10)  were  used  in  the  calculation  of  a planar 
diffuser  flow  field  with  a fully  developed  parabolic  entrance  velocity 
profile.  The  detailed  numerical  results  for  separated  and  non-separated 
flows,  along  with  an  evaluation  of  convergence  characteristics  of  the  nu- 
merical method  used,  are  discussed  in  Section  6. 1. 


32  ALGEBRAIC  VISCOSITY  MODEL 

In  the  algebraic  viscosity  model  for  boundary -layer  flows,  the  eddy 
viscosity  is  related  to  the  mean  velocity  field  parameters  through  an 
algebraic  relationship.  There  are  two  types  of  algebraic  models  com- 
monly used,  namely,  a local  model  and  a global  model.  The  former, 
such  as  Prandtl's  mixing  length  theory,  relates  the  eddy  viscosity  to 
the  local  velocity  gradient.  On  the  other  hand,  a global  model  such  as 
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the  one  proposed  by  Clauser  (Eef.  15)  relates  the  eddy  viscosity  to  the 
integral  quantities  such  as  the  displacement  thickness  of  a boundary 
layer.  The  Prandtl's  mixing  length  model  is 


(16) 


where  S.  is  the  mixing  length.  In  the  wall  region,  the  mixing  length  is 
given  by  Van  Driest  (Ref.  16)  as 

*=0.  *l»[i-exp(-j*7Jfr)]  (17) 

where  tw  is  the  shear  stress  at  the  wall. 

The  Van  Driest  model  has  been  widely  used  and  modified  to  include 
other  effects  such  as  the  pressure  gradient  (9p/9x),  wall  suction,  etc. 
(Ref.  17).  The  model  has  been  used  to  obtain  fairly  good  results  for 
certain  boundary-layer  flows  (Refs.  18  through  21).  It  is  used  in  the 
present  analysis  to  calculate  the  velocity  distribution  for  a fully  devel- 
oped channel  flow.  The  fully  developed  velocity  profile  is  then  used  as 
the  inlet  condition  for  a diffuser  flow  calculation.  For  fully  developed 
channel  flows,  the  total  shear  stress  can  be  written  as 


(18) 


With  the  assumption  that  shear  stress  is  constant  near  the  wall, 
i.e.,  t *tw>  Eqs.  (16),  (17),  and  (18)  are  used  to  derive  the  velocity 
gradient: 


where 


2 

I +J I + 4 c 0.40*  y*1  [ I - exp  (-  2 


(19) 


and 


(20) 
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From  the  vorticity  distribution,  given  by  Eq.  (19),  the  velocity  and  the 
stream  function  distributions  can  be  easily  obtained  by  numerical  inte- 
gration. 

For  more  complicated  turbulent  flows,  such  as  those  with  strong 
adverse  pressure  gradients  and  with  separation,  Prandtl's  mixing 
length  model  is  not  adequate  to  obtain  satisfactory  results.  Although 
extensions  and  modifications  of  the  model  are  possible,  they  usually  re- 
quire experience  in  the  turbulence  modeling  of  a particular  flow  problem. 
Often,  a trial-and-error  approach  or  a numerical  optimization  proce- 
dure is  employed  to  obtain  the  desired  eddy  viscosity  distribution. 


3.3  DIFFERENTIAL  EQUATION  VISCOSITY  MODEL 

In  the  differential  equation  models,  the  eddy  viscosity  (i^)  is  re- 
lated to  the  characteristics  of  the  turbulent  motion,  such  as  the  turbu- 
lent kinetic  energy,  which  are  obtained  by  solution  of  additional  differ- 
ential equations.  Two  models  in  this  category  are  commonly  used, 
namely,  (1)  the  Prandtl's  one -equation  turbulent  kinetic  energy  model 
and  (2)  the  Prandtl -Kolmogorov  two-equation  model  (Ref.  13). 

For  the  Prandtl's  turbulent  kinetic  energy  model,  the  eddy  viscosity 
is  related  to  the  product  of  a length  scale  and  the  square  root  of  the  tur- 
bulent kinetic  energy,  i.  e.  , 


^ B l’fi<  (21) 

where  H is  a length  scale.  The  turbulent  kinetic  energy  (k)  is  obtained 
from  a differential  equation  and  the  length  scale  ($.)  is  specified  alge- 
braically. For  simple  turbulent  flows,  several  length  scale  formula 
have  been  proposed  (Ref.  13).  Unfortunately,  the  specification  of  the 
length  scale  still  requires  experience  and,  more  often,  a trial-and- 
error  approach.  Therefore,  the  method  is  not  general  and  cannot  be 
used  for  complex  turbulent  flows  with  ease.  Nevertheless,  Prandtl's 
turbulent  kinetic  energy  model  contains  more  information  about  the 
turbulent  motion  than  does  his  mixing  length  theory. 

The  TKE  equation,  which  can  be  derived  from  the  Navier -Stokes 
equation  (Ref.  12),  is  written  as 
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u.jk  ..JLf.y 
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- U^U' 

L J 


_ / ?u[  ^U[  > 

0Xj  ~ V'2Xj  dXj  ' 


where  the  summation  convention  has  been  used. 

Equation  (22)  contains  convection,  diffusion,  and  source  terms. 
Because  the  turbulent  kinetic  energy  is  dependent  on  the  history  of  the 
turbulence  motion,  the  eddy  viscosity  model  (Eq.  (21))  also  depends  on 
the  history  of  the  turbulence  motion.  Examples  of  turbulent  flows  with 
a strong  history  dependent  nature  are  flows  with  separations,  flows 
with  strong  adverse  pressure  gradient,  etc. 

For  the  Prandtl -Kolmogorov  two-equation  model,  the  length  scale 
(j2)  is  obtained  from  a transport  equation.  For  high  Reynolds  number 
regions  of  the  flow  where  the  molecular  viscosity  effect  is  small,  the 
length  scale  associated  with  the  turbulent  kinetic  energy  dissipation 
(e  = (k)3/2/I  e)  is  used  in  Eq.  (21)  to  obtain 


**  A 8 


(23) 


where  is  assumed  to  be  a constant  (0.09)  at  high  Reynolds  number. 
Although  additional  equations  for  e and  k are  needed,  the  specification 
of  the  length  scale  is  completely  eliminated  which  makes  the  model  at- 
tractive for  the  calculation  of  complex  turbulent  flows. 


3.3.1  High  Reynolds  Number  Two-Equation  k-e  Model 

For  high  Reynolds  number,  the  governing  equations  for  the  turbu- 
lent kinetic  energy  and  its  dissipation  are  modeled  as  (Ref.  22) 

k- Equation 


(24) 
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e -Equation 


+ 

where 


(25) 


c*,  = J.44  , = 1.92  and  &£  =,.| 

The  value  of  constants  used  follows  that  of  Launder  (Ref.  22),  and  e is 
the  so-called  isotropic  part  of  the  total  TKE  dissipation. 


In  deriving  Eqs.  (24)  and  (25),  it  was  assumed  that  the  effect  of  the 
molecular  viscosity  is  negligible.  The  assumption  is  valid  in  the  fully 
turbulent  regions  at  high  Reynolds  number  where  the  turbulent  eddy  vis- 
cosity (vj)  is  usually  several  order  of  magnitude  larger  than  the  molec- 
ular viscosity  ( v ).  Inside  the  viscous  wall  sublayer,  however,  the  as- 
sumption is  no  longer  valid  because  the  molecular  viscosity  does  play  a 
dominant  role  there.  For  this  reason,  a special  sublayer  formulation 
is  necessary  which  can  be  accomplished,  for  example,  by  matching  the 
solution  from  the  fully  turbulent  region  to  an  analytical  law  of  the  wall 
solution.  Two  types  of  analytical  expressions  are  commonly  used, 
namely,  a logarithmic  law  of  the  wall  and  a power  law  velocity  profile. 
The  former  is  given  as 

U = V#[2.5(«(^VW+5.5]  <26> 

where  y is  the  distance  from  the  wall  and  v*  is  the  friction  velocity  de- 
fined as  \F^7p. 


In  the  present  study,  the  high  Reynolds  number  k-e  model  with  the 
law  of  the  wall  matching  is  used  in  the  calculation  of  a conical  diffuser 
flow  field  with  a fully  developed  turbulent  pipe  flow  at  the  entrance. 

The  results  are  presented  in  Section  6.3. 


3.3.2  Low  Reynolds  Number  Two-Equation  k-e  Model 

For  turbulent  separated  diffuser  flow  calculations,  the  prediction 
of  the  point  of  separation  and  the  reverse  flow  field  in  the  separated  re- 
gion is  of  major  importance.  The  high  Reynolds  number  k-e  model 
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with  the  law  of  the  wall  matching  procedure  is  not  applicable  in  the  sep- 
arated flow  calculation  because  Eq.  (26)  becomes  singular  at  the  point 

of  separation,  i.  e. , v*  = \ rw/ p = 0,  In  addition  the  validity  of  the  law 
of  the  wall  profile  (Eq.  (26))  in  a flow  with  a strong  adverse  pressure 
gradient  is  also  questionable  (Refs.  23  and  24).  Therefore,  it  is  neces- 
sary to  include  the  sublayer  in  the  formulation  so  that  the  point  of  sep- 
aration can  be  predicted.  In  the  present  analysis,  this  is  accomplished 
by  the  use  of  (1)  a low  Reynolds  number  two-equation  k-e  model  that 
includes  the  effect  of  the  molecular  viscosity,  and  (2)  a coordinate 
transformation  which  stretches  the  sublayer  region  to  provide  many 
computational  grid  points  in  that  region. 

For  low  Reynolds  number,  in  Eq.  (23)  is  redefined  in  the  pres- 
ent study  as 

" Ica-hA/h)  (27) 


where  A = y[2ky !v,  y is  measured  from  the  wall,  and  a and  b are  con- 
stants. The  k- equation  becomes 


/*♦  *1- -1—  {rtt-  £Sh  A +rv- 24. 

tsx*  9r‘J  0>+>plL  '3xJ9x^Lv  it  r hrf 


(28) 


The  last  term  of  Eq.  (28)  is  the  total  TKE  dissipation  which  consists  of 
the  isotropic  and  the  non -isotropic  parts  of  the  dissipation.  The 
e- equation  becomes 


l«x*  ar‘i  ?x  +iv  ar'<r£)  r 


(29) 


I 


o 
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where 

C?,  = 1.44  , = f.32£i-  o.iexp(-R1)]  , -|.|  (30) 

and' 

R=  k*/(j>£ ) (3D 

Note  that  Eqs.  (27),  (28),  and  (29)  reduce  to  Eqs.  (23),  (24),  and  (25), 
respectively,  when  v is  neglected.  The  low  Reynolds  number  terms  in 
Eqs.  (27),  (28),  and  (29)  are  derived  in  the  following  paragraphs. 

3.3.2.1  Derivation  of  the  Total  Turbulent  Kinetic 
Energy  Dissipation  at  the  Wall 


The  exact  expression  for  the  total  turbulent  kinetic  energy  dissipa- 
tion is  (Ref.  12): 


<r  - & 'ft ♦#♦<#♦#*#  ♦£/}  <»> 


At  the  wall,  velocity  gradients  such  as  9u'/9x,  9u'/9z,  9v'/9x,  9v'/9z, 
9w'/9x,  and  9w'/9z  vanish.  In  addition,  from  the  continuity  equation, 
9v'/9y  also  vanishes.  Therefore,  the  total  turbulent  kinetic  energy 
dissipation  at  the  wall  can  be  simplified  to 


(33) 


Equation  (33)  can  be  evaluated  at  the  wall  in  the  following  manner: 


(34) 


where  y is  measured  from  the  wall.  Equation  (34)  can  also  be  written 
in  terms  of  the  TKE  as 
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<3! 

The  last  term  in  Eq.  {35)  is  small  compared  to  the  first  term  near  the 
wall  and,  therefore,  may  be  neglected.  The  final  expression  for  the 
total  turbulent  kinetic  energy  dissipation  at  the  wall  becomes 


Z V 


<36) 


This  term  was  ignored  in  the  high  Reynolds  number  formulation  of  the 
turbulent  kinetic  energy  equation  (Eq.  (24)).  But,  it  should  be  included 
in  the  low  Reynolds  number  formulation  in  regions  where  the  molecular 
viscosity  effect  is  not  negligible.  In  the  high  Reynolds  number  regions, 
the  total  TKE  dissipation  consists  of  only  the  isotropic  part  of  the  dis- 
sipation. The  isotropic  part  decreases  rapidly  near  the  wall  and  van- 
ishes completely  at  the  wall.  Since  the  total  TKE  dissipation  term  is 
needed  in  the  calculation  of  the  low  Reynolds  number  TKE  equation  not 
only  at  the  wall  but  also  throughout  the  viscous  sublayer,  a formula  that 
represents  the  non -isotropic  part  of  the  turbulent  kinetic  energy  dissi- 
pation needs  to  be  developed.  In  the  present  study,  Eq.  (3  6)  is  used  not 
only  at  the  wall  to  represent  the  TKE  dissipation  but  also  throughout  the 
flow  field  to  represent  the  non-isotropic  part  of  the  TKE  dissipation. 

The  total  TKE  dissipation  is  represented  by 

£t  <37> 

throughout  the  flow  field.  Equation  (37)  appears  as  the  last  term  in  the 
low  Reynolds  number  formulation  (Eq.  (28)). 

3.3.2.2  Derivation  of  the  Coefficient  (C^)  for  the 
Low  Reynolds  Number  Model 

In  order  to  determine  the  coefficient  (C^)  from  Eq.  (2*3),  one  needs 
to  have  a realistic  distribution  of  and  k2/e  near  the  wall.  Since  one 
of  the  ultimate  goals  of  the  present  analysis  is  the  prediction  of  sep- 
arated diffuser  flows,  the  eddy  viscosity  models  proposed  by  Mellon 
and  Herring  (Ref . 25),  Glusko  (Ref.  26),  and  Alber  (Ref.  27)  for  nearly 
separated  or  separated  boundary -layer  flows  are  used.  These  authors 
represent  the  eddy  viscosity  in  the  sublayer  region  by 
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*t  - v( 


CL+  (A/iO 


(38) 


where  A is  the  characteristic  length  in  the  sublayer  defined  as  y}  2k  y/v 
and  a and  b are  constants.  The  total  turbulent  kinetic  energy  dissipa- 
tion is  modeled  by  Mellon  and  Herring  as 


9 A2 

VTF  (i+A/ti) 


I 


(39) 


The  isotropic  part  of  the  TKE  dissipation  can  be  obtained  from  Eqs. 
(39)  and  (37)  as 

£ = k */(}>>  A) 

Equation  (40)  can  also  be  written  as 


(41) 


By  substituting  Eqs.  (41)  and  (38)  into  Eq.  (23),  the  final  expression  for 
the  coefficient  (C^)  is 


£ - 0.  - 
s1  3(a+A/b) 


(42) 


In  this  study,  the  value  of  a is  taken  to  be  1,100,  while  the  value  of  b is 
determined  by  taking  the  high  Reynolds  number  limit  of  C^,  i.  e.  , 


5“ 


A -*  Large 


y = 009 


(43) 


where  the  value  0. 09  appears  to  be  representative  for  high  Reynolds 
number  flows  (Refs.  28  and  29).  Thus,  the  coefficient  (b)  determined 
from  Eq.  (43)  is 

b a 0.27  (44) 

The  function  is  plotted  in  Fig.  1.  It  should  be  noted  that  direct  mea- 
surement of  the  dissipation  term  e is  difficult.  Thus,  the  validity  of  the 
analytical  model  expressions  can  be  verified  only  through  the  comparison 
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of  the  numerically  computed  and  experimentally  measurable  quantities 
in  turbulent  flows. 

3.3.3  Discussions  on  the  Present  and  Launders  Model 


There  are  some  basic  differences  between  the  present  and 
Launder's  model.  The  approach  used  by  Launder,  et  al. , in  determin- 
ing the  coefficient  (C^)  was  based  on  the  less  general  Van  Driest  form 
of  the  mixing  length  formula  (Eq.  (17))  which  determined  the  eddy  vis- 
cosity distribution  in  a constant  shear  Couette  flow.  The  eddy  viscosity, 
thus  determined,  was  used  on  the  left-hand  side  of  Eq.  (23).  The 
e -equation  was  then  "adjusted  and  calculated"  so  that  a reasonable  tur- 
bulent kinetic  energy  distribution  was  obtained  in  the  viscous  sublayer 
region.  Finally,  with  k and  e determined,  Eq.  (23)  was  inverted  to 
provide  a preliminary  estimate  of  Cju.  The  functional  form  of  C ^ ob- 
tained from  Eq.  (23)  was  then  used  to  provide  the  final  expression  which 
is  given  in  Table  1 . 


During  the  course  of  the  numerical  optimization,  Launder  found 
that  an  additional  artificial  term  is  needed  in  the  e -equation  in  order  to 
have  a maximum  value  of  the  TKE  distribution  at  y+  « 20  as  indicated 
by  experimental  data.  The  final  form  of  Launder's  artificial  term  is 
also  given  in  Table  1 . 

Consider  the  modeling  of  the  non-isotropic  part  of  TKE  dissipation 
term.  From  Table  1,  it  can  be  seen  that  both  models  use  the  TKE  dis- 
tribution (k).  For  the  purpose  of  demonstration,  a typical  TKE  distri- 
bution in  a fully  developed  channel  flow  is  given  in  Fig.  2.  The  turbu- 
lent kinetic  energy  reaches  a maximum  value  near  the  wall  and  decreases 
monotonically  toward  the  wall  and  the  centerline.  Since  the  present  dis- 
sipation model  is  directly  proportional  to  the  turbulent  kinetic  energy, 
it  is  always  positive.  The  present  model  provides  a non-zero  value  at 
the  maximum  TKE  location  (yj"  ).  Because  the  slope  of  the  TKE  dis- 

tribution  vanishes  at  the  peak  location,  the  Launder's  model  for  (e^  - e) 
is  zero  (see  Fig.  2).  As  a result,  Launder's  model  underestimates  the 
energy  dissipation  in  the  neighborhood  of  y^ax  which  probably  affects  the 
solution  in  such  a way  that  both  the  TKE  and  its  second  derivative  can- 
not reach  the  experimental  value. 


With  the  present  model,  Eq.  (28)  becomes 


9rl  (?+>>)  ] 


(V 


7^{£+"'-F}-0 


(45) 


24 


AE  DC-T  R-76-1 69 


With  the  Launder's  model,  Eq,  (28)  becomes 


ik 

9F1 


Pt  (4  Uf  I 


= 0 


(46) 


The  last  term  of  Eq.  (46)  represents  only  the  isotropic  part  of  the 
energy  dissipation.  To  overcome  this  difficulty.  Launder  created  an 
artificial  term  (Table  1)  in  the  e -equation  so  that  the  magnitude  of  the 
isotropic  part  of  TKE  dissipation  (e)  can  be  increased  to  a higher  level. 


4.0  COORDINATE  SYSTEMS  AND  TRANSFORMATIONS 


The  governing  equations  introduced  in  the  previous  sections  must 
be  solved  in  conjunction  with  boundary  conditions  for  the  particular  flow 
problem  of  interest.  The  key  to  the  successful  use  of  a given  formula- 
tion is  strongly  a function  of  the  coordinate  system  used  and  the  nu- 
merical method  adopted  in  the  calculation.  Several  types  of  coordinate 
systems  and  transformations  may  be  used  depending  on  the  problem  to 
be  solved.  Since  the  governing  equations  and  boundary  conditions 
usually  are  discretized  in  finite  difference  form  and  then  solved  alge- 
braically on  digital  computers,  important  factors,  such  as  the  computer 
storage,  the  computing  time,  and  the  accuracy  of  the  solution,  need  to 
be  considered  in  the  formulation  of  a computer  program.  The  size  of 
the  computer  storage  limits  the  number  of  discretized  variables  one 
can  use  for  a particular  flow  problem.  It  also  limits  the  accuracy  of 
the  solution  because  of  the  finite  number  of  grid  points  available  to  de- 
scribe the  flow  field.  Naturally,  the  more  the  grid  points  the  higher  the 
accuracy. 


4.1  UNIFORM  AND  NONUNIFORM  SYSTEMS 

The  simplest  way  to  describe  a flow  field  is  to  use  a uniform  co- 
ordinate system  that  has  equal  spacing  in  each  of  its  coordinates  (see 
Fig.  3).  The  advantages  are  (1)  the  physical  location  can  be  easily 
identified  which  eases  the  interpretation  of  the  solution,  and  (2)  fewer 
calculations  are  needed  because  the  grid  spacing  needs  to  be  calculated 
only  once  in  a program.  In  the  present  study,  the  uniform  system  is 
used  in  Section  6. 1 to  calculate  the  flow  field  in  a diffuser  with  a para- 
bolic inlet  velocity  profile.  The  diffuser  geometry  used  has  constant 
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width  inlet  and  exit  sections  which  fit  the  coordinate  system  nicely  ex- 
cept at  the  boundary  of  the  diverging  section.  The  diverging  section  can 
be  represented  by  generating  the  coordinate  lines  from  the  diverging 
section  so  that  grid  points  lie  along  the  diverging  wall  (see  Fig.  3b). 

This  approach  works  well  when  the  diverging  angle  is  neither  too  large 
nor  too  small. 

When  large  velocity  gradients  occur  in  the  flow  field,  the  uniform 
system  is  no  longer  adequate  to  describe  the  flow.  In  that  situation,  a 
nonuniform  system  is  necessary.  By  arranging  more  grid  lines  in  the 
region  where  large  velocity  gradients  occur  (Fig.  3c),  the  accuracy  of 
the  solution  can  be  improved.  But  the  approach  requires  experience  on 
a particular  flow  problem,  and  in  most  cases,  the  manual  arrangement 
of  the  coordinates  is  inevitable.  In  addition,  the  accuracy  of  the  re- 
sult from  a nonuniform  system  is  difficult  to  interpret  in  terms  of  grid 
spacings.  Nevertheless,  the  nonuniform  coordinate  system  is  highly 
flexible  in  the  early  stages  of  program  development  so  that  the  basic 
features  of  the  solution  can  be  obtained.  This  approach  was  adopted 
for  the  calculation  of  a turbulent  separated  diffuser  flow  with  an  algebraic 
turbulence  model  in  Section  6.2. 


4.2  COORDINATE  TRANSFORMATION 
4.2.1  Body-Aligned  Coordinate  Transformation 


For  diffusers  with  a small  diverging  angle,  it  becomes  very  diffi- 
cult to  use  the  nonuniform  coordinate  system  such  as  the  one  shown  in 
Fig.  3c.  A body-aligned  coordinate  system  is  essential  for  a proper 
solution  in  this  case,  and  it  also  works  well  for  curved  walls.  The  co- 
ordinate transformation  used  in  the  present  analysis  is 


X 


* X 


(47) 


where  S(x)  represents  the  lateral  coordinate  of  the  diffuser  wall.  Equa- 
tion (47)  is  a linear  stretching  function.  The  corresponding  transforma- 
tion factors  are 


9x  _ . 9x  _ 0 2L - J ar  _ - £<:*) 

?x  * ar  ” ' ar  "S(x)  * ax  “”r  g#) 


(48) 


where  S'(x)  represent  dS(x)/dx. 


The  transformation  given  by  Eq.  (47) 


26 


AEDC-TR-76-159 


transforms  a physical  domain  of  a diffuser  into  a rectangular  shape  (see 
Fig.  4).  The  governing  equations  developed  in  the  previous  sections  are 
transformed  into  the  x,  r plane  by  using  the  chain  rule  differentiation. 


m 9 (3L)1 

flx  9x 


JL  - /JL\  JL 
ir  " { -ar  ;9r 


9X  _ V 


+ ( 


9 r 

9X 


)*— 

;,p» 


9*  _ 9 *L 

9xar"9xar  (»r;i3x^  9?1 


jl./iLfjL 

> ar2 


(49) 


Although  the  transformed  domain  shape  becomes  simple,  the  num- 
ber of  terms  in  the  governing  equations  grows.  In  the  present  analysis, 
Eq.  (47)  is  used  to  transform  an  8-deg  conical  diffuser  into  a rectangu- 
lar shape.  The  computation  was  then  carried  out  in  the  transformed 
domain.  The  results  are  presented  in  Section  6.3. 

4.2.2  Coordinate  Transformation  with  a Sublayer  Stretching 


The  existence  of  a thin  sublayer  commonly  found  in  the  turbulent 
wall  boundary  layer  makes  it  very  difficult  to  solve  the  whole  problem 
in  the  physical  domain  with  a uniform  coordinate  system.  It  becomes 
necessary  to  stretch  the  sublayer  region  in  such  a way  that  the  sharp 
gradient  in  flow  variables  can  be  reasonably  resolved.  In  the  present 
analysis,  a composite  coordinate  transformation  is  used  to  provide 
good  resolution  throughout  the  flow  field.  For  a fully  developed  channel 
flow  velocity  profile,  the  sublayer  velocity  distribution  is  a linear  func- 
tion of  the  distance  from  the  wall,  i.  e. , 


1 L*  = / (50) 

where  u+  = u/v*  and  y+  = y v*/v.  The  velocity  gradient  of  the  sublayer 
profile,  i.  e.  , 9u/9y  = v'^/v,  provides  the  basic  stretching  factor  for 
9y / 9y  in  the  coordinate  transformation  given  by 


U = J P?  , 0 « S S a. 


(51) 


where  y is  the  coordinate  measured  from  the  wall  with  y = 1 at  the 
centerline,  y is  the  stretched  coordinate,  and  a,  j3,  and  yQ  are  parameters 
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in  the  transformation.  In  the  core  region,  the  stretching  factor  (ay/ay) 
is  gradually  changed  to  unity  at  the  centerline  where  no  stretching  is 
needed  because  of  the  low  gradients  in  the  velocity  profile.  The  smooth 
transition  is  provided  by  the  core  transformation  function 

<f  = Cjf  + Irvfcosl,  (»-»„)]  +F  , (52) 

where  c and  F are  parameters  in  the  transformation.  By  proper  match- 
ing of  the  two  transformation  functions,  one  can  determine  the  coeffi- 
cients and  thus  provide  a continuous  transformation  for  the  whole  flow 
field.  A typical  coordinate  transformation  given  by  Eqs.  (51)  and  (52) 
is  shown  in  Fig.  5.  The  detailed  derivation  of  the  parameters  used  in 
the  transformation  is  given  in  Appendix  A. 

The  ability  of  the  coordinate  transformation  to  provide  good  resolu- 
tion in  the  transformed  coordinate  plane  is  demonstrated  by  considering 
the  velocity  and  the  turbulent  kinetic  energy  distributions  in  a channel 
flow  shown  in  physical  coordinates  in  Fig.  6a.  It  can  be  seen  that  large 
gradients  of  the  velocity  and  turbulent  kinetic  energy  profiles  exist  near 
the  wall.  Poor  resolution  can  be  expected  in  the  sublayer  region  when 
one  attempts  to  use  a uniform  coordinate  system  to  describe  the  pro- 
files. However,  when  Eqs.  (51)  and  (52)  are  used  to  stretch  the  sub- 
layer and  the  core  region,  sharp  gradients  in  velocity  profiles  diminish 
as  indicated  in  Fig.  6b.  Thus,  the  detail  velocity  and  turbulent  kinetic 
energy  profiles  can  be  adequately  described  in  the  transformed  coordi- 
nates. For  this  reason,  the  coordinate  transformation  plays  an  impor- 
tant role  in  obtaining  an  accurate  numerical  solution. 

4.2.3  A Complete  Coordinate  Transformation  for  a Diffuser  Flow 

Both  the  body-aligned  coordinate  transformation  and  the  sublayer 
stretching  are  necessary  to  provide  a good  spatial  resolution  through- 
out a diffuser  flow  field.  The  complete  coordinate  transformation  is 
achieved  first  by  mapping  the  diffuser  shape  into  a rectangular  domain 
then  followed  by  a sublayer  stretching  (see  Fig.  7).  The  complete 
transformation  is  given  as 

X=X  (53) 


In  the  sublayer  region. 


*"ma* 


(54) 
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In  the  core  region, 

r=  $ («>{i  - [c  c r ) f £„[  cosh  + F ]} . 


<55) 


_ - 

osr^r, 


where  r is  measured  from  the  centerline,  S(x)  represents  the  diffuser 
wall  shape,  rQ  is  the  transformed  matching  location,  and  rmax  is  the 
transformed  wall  shape.  The  complete  transformation  factors  are: 

In  the  sublayer  region  from  Eq.  (54), 
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In  the  core  region,  from  Eq.  (55), 
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(57) 
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The  transformation  factors  given  by  Eqs.  (56)  and  (57)  will  be  used  in 
the  next  section  to  derive  the  governing  equations  in  the  transformed 
coordinates.  In  the  present  analysis,  the  coordinate  transformation 
given  by  Eqs.  (53),  (54),  and  (55)  was  used  in  the  calculation  of  sep- 
arated and  non-separated  diffuser  flows.  The  results  are  presented 
in  Section  6.4. 


5.0  FINITE  DIFFERENCE  FORMULATION  AND  NUMERICAL 
SOLUTION  PROCEDURE 


The  governing  equations  and  the  turbulence  models  presented  in 
the  previous  sections  are  a system  of  coupled  non-linear  partial  differ- 
ential equations.  The  system  cannot  be  solved  analytically  and,  there- 
fore, must  be  solved  by  numerical  methods.  The  system  of  equations 
becomes  even  more  complicated  when  it  is  written  in  the  transformed 
coordinates  such  as  those  described  in  Section  4.  2.  3.  In  the  present 
analysis,  a standard  form  of  the  transformed  equations  is  first  de- 
rived to  represent  the  common  features  of  the  governing  equations.  A 
general  finite  difference  formulation  is  then  developed  so  that  stable 
and  convergent  solutions  can  be  obtained  for  a wide  range  of  Reynolds 
number.  The  stability  limitation  associated  with  the  central  difference 
scheme  and  the  accuracy  problem  inherent  to  the  up-wind  difference 
scheme  are  avoided  by  the  use  of  the  locally  evaluated  decay  functions 
in  the  finite -difference  formulation. 


5.1  COMPLETE  GOVERNING  EQUATIONS  IN  THE  TRANSFORMED 
COORDINATES  AND  THE  STANDARD  FORM  EQUATION 

For  the  separated  or  non-separated  diffuser  flow  calculations, 
the  coordinate  transformation  given  by  Eqs.  (53),  (54),  and  (55)  is 
necessary  to  provide  adequate  resolution  in  both  the  core  and  the  vis- 
cous sublayer  region.  The  vorticity-stream  function  formulation  and 
the  two-equation  low  Reynolds  number  turbulence  model  written  in  the 
transformed  coordinates  are: 
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Stream  Function  (y;): 

± - l£r- ££  tL 1 & J \ . , / jL  j^C 

9X1-  ariax^'arM  5?  7 * Sr ' + z * a* ; aSfop 


. A»*  T _ 

9xl  «r 

Velocity  Recovery  (u,  v): 


P + r a 


U=  -L(i£)2£ 

u r«  'jr'a? 

Lr23e+5Sf(2£)l 

Turbulent  Kinetic  Energy  Equation  (k): 

..  

- {tv-  £ (£)  4 < **>1  fg>  ♦ [u-  (£+ 

+ f]‘+r^#]*+*^l}+[?F^)+,T+ <&#rf 

k“° 

where  y is  measured  from  the  wall. 


0 

(59) 

(60) 
(61) 


(62) 
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Turbulent  Kinetic  Energy  Dissipation  (e): 


<s / I 


• (y+ W l *2  t 

^ i-atf  /»l\  * 2L  wl£}2112l  - d 5 « 

+ [fir  (9r'+W+(**J9r  J J 1 k (?+*/& 


(63) 


The  Prandtl-Kolmogorov  eddy  viscosity  formula  is 


The  coefficients  (C^,  C2,  cxe,  and  C^)  are  given  as 

£,=  1.44  , Ct  = f.92Ci-o.3exp(-iRa)  ] , «£  = |.l 
- A/C3  (1100  +A/0.2  7)] 

where 


(64) 


(65) 


The  equations  are  lengthy  and  complicated,  especially,  when  the 
transformation  factors  are  calculated  from  Eqs.  (56)  and  (57).  To  sim- 
plify the  analysis,  Eqs.  (58),  (59),  (62),  and  (63)  can  be  cast  in  a 
standard  form  which  retains  the  basic  features  such  as  the  diffusion, 
the  convection,  the  production  and  the  dissipation  of  a flow  variable. 

The  standard  form  derived  is 


where  <j>  represents  the  flow  variables,  i.  e. , y/,  k,  and  e.  The  cor- 

responding coefficients  (a^,  a.2,  bj,  b2,  and  d)  are  given  below  for  each 
of  the  variables,  , 1 p , k,  and  e.  ' 


32 


AE  DC-T  R-76-1 59 


-equation: 

Q,=  I 
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k-equation: 

o,=  l 
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e - equation: 
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5.2  A GENERAL  FINITE  DIFFERENCE  FORMULATION  WITH 
DECAY  FUNCTIONS 


The  partial  differential  equation  (Eq.  (66))  can  be  written  as 

L{4>}  +d  = 0 (71) 

where  L represents  a differential  operator.  Equation  (71)  written  in 
a conventional  finite  difference  form  is 


Uj.{4>}  = e (72) 

where  is  a finite -difference  operator  and  e represents  the 
difference 


In  the  present  approach  (Ref.  30),  the  finite -difference  equation  is 
written  as 


+<1f  a » 


(74) 


where  the  additional  function  (G)  is  named  as  the  "decay  function" 
which  is  used  to  minimize  the  truncation  error  (7).  The  expression 
for  the  decay  function  (G)  is  derived  in  Appendix  B. 


The  present  finite  difference  formulation  of  the  standard  Eq. 
(66)  is 


a f ^i>l > j " 2 > J + J I / b|  x } 

a,l IP eTW  *»*  I 

* + ai.j  -0 


where  the  decay  functions  (Gj  and  Gj)  are  determined  from 

^t»(-|:)ri-2(eRt-i)]  (eRl-0 


(75) 


(76) 
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< 77) 


and  the  grid  Reynolds  numbers  (R^  and  Rj)  are  defined  as 


(78) 


The  decay  function  is  shown  in  Fig.  8.  For  the  conventional  central 
difference  scheme,  the  decay  functions  (Gi  and  Gj)  are  set  equal  to 
unity,  i.  e. , 


„ f Vr2*i,j4<W  / b,  % 4ur4Uj  i 
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32): 


The  well-known  stability  limitation  for  Eq.  (79)  is  (Refs.  31  and 


sr  £2 


(80) 


which  for  the  conventional  central  difference  scheme  severely  limits 
the  grid  sizes  (6x  and  6r)  when  the  coefficients  (bj/a^  and  b2/a2)  are 
very  large.  In  order  to  obtain  a stable  and  convergence  solution,  the 
total  number  of  grid  must  be  increased  when  the  small  grid  meshes 
are  needed,  which  creates  computer  core  storage  problems.  The 
stability  problem  is  avoided  by  using  the  decay  functions  given  by 
Eqs.  (76)  and  (77).  Equation  (75)  is  written  as 


, A|  X f Vi  ~ZK'i  + Vi  / b.  /■»  \ Vi  ■ Vj  1 
{Qi][  s*2  HX  J 

+ / f V+'  “2Vi  + Vrl  _ / bi  G % l + J . . - A 

4 VI TF J + d^“0 


(81) 


Equation  (81)  is  similar  to  Eq.  (79)  except  that  the  coefficient 
(a!  is  replaced  by  a.1/Gi  and  bj/a!  is  replaced  by  (bi/ai)Gi,  etc. 
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Therefore,  the  stability  limitation  (Eq.  (80))  can  be  written  as 
After  re-arrangement,  Eq.  (82)  becomes 


> 


(83) 


which  is  the  upper  bound  for  and  Gj  shown  in  Fig.  8.  It  can  be  seen 
from  Fig.  8 that  the  decay  functions  (Gj)  given  by  Eq.  (77)  is  lower 
than  the  upper  bound  of  the  stability  limitation.  Therefore,  there  is 
no  stability  limitation  on  the  grid  size  (6x  or  &T)  with  the  present 
general  finite  difference  formulation  because  Eq.  (83)  is  automatically 
satisfied  with  G^,  Gj  determined  from  Eqs.  (76)  and  (77).  The  uncon- 
ditional stability  is  achieved  because  additional  analytical  information 
has  been  incorporated  into  the  derivation  of  decay  functions  (Gi  and 
Gj)  (Appendix  B).  This  characteristic  of  the  present  formulation  makes 
it  possible  to  solve  the  complicated  equations  given  by  Eqs.  (66) 
through  (70)  and  obtain  stable  convergence  solutions. 

The  calculation  of  the  decay  functions  (Gi  and  Gj)  can  be  simplified 
by  the  following  approximation. 

G = I.O-  0,625  „ IRK  2 

= -1 L_  m]>z  <84> 

IRI  CR)1  > 

for  both  Gj  and  Gj,  and  Rj  and  Rj.  The  approximation  (Eq.  (84))  is 
also  shown  in  Fig.  8. 


It  is  also  interesting  to  examine  the  high  grid  Reynolds  number 
limit  of  the  present  formulation.  When  the  grid  Reynolds  number  is 
large,  the  decay  function  can  be  approximated  as 


(85) 
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By  using  Eq.  (85),  the  finite  difference  formulation  (Eq.  (81))  becomes 


> l^i.1  . f / bi  2 \ J j 1 

1 2.  Bx * ' K ai  IRj.1  *6*  I 

/lRila»x  f ^j+» - 2^,j  tb,  1 \ ? . J 

+ ( 2 M sp (a*  Tr-t  IS? — J + d^js0 


(86) 


For  Rp  Rj  > 0 and  by  using  Eq.  (78),  Eq.  (86)  can  be  written  as 


(87) 


Equation  (87)  is  the  up-wind  difference  formulation  of  the  high  grid 
Reynolds  number  inviscid  equation: 

'{b'  fr  + b*sf}  + d *°  <8B) 

which  may  be  derived  directly  from  Eq.  (66)  by  neglecting  the  viscous 
terms.  Note  that  Eq.  (88)  is  valid  outside  the  boundary  layer  region. 
Thus,  the  present  finite  difference  formulation  approaches  the  con- 
ventional central  difference  scheme  at  low  grid  Reynolds  number  but 
gradually  changes  into  the  up-wind  difference  scheme  as  the  local 
grid  Reynolds  number  (R^Rj)  increases.  This  transition  is  provided 
by  the  decay  functions  (G^  and  Gj),  which  are  determined  through  the 
use  of  local  analytical  solutions  (Appendix  B). 


5.3  NUMERICAL  SOLUTION  PROCEDURE 

When  the  finite  difference  equation  (Eq.  (75))  is  applied  at  each 
grid  point  for  each  flow  variable,  a large  system  of  algebraic  equa- 
tions is  formed.  The  system  of  equations  and  the  corresponding 
boundary  conditions  are  solved  in  the  present  analysis  by  a standard 
Gauss-Seidel  point -iteration  scheme  (Refs.  33  and  34).  This  approach 
has  several  advantages,  namely,  (1)  a relatively  small  amount  of 
computer  storage  is  required  (one  location  for  each  variable),  (2)  the 
program  is  easily  written,  and  (3)  under-  or  over -relaxations  can  be 
easily  incorporated  so  that  the  rate  of  convergence  can  be  optimized. 
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The  general  solution  procedure  is  shown  in  Fig.  9.  The  iteration  pro- 
cedure begins  with  an  initial  guess  of  the  flow  field  which  can  be  obtained 
by  assuming  that  the  flow  profiles  are  similar  to  the  inlet  profiles.  The 
inlet  profiles  must  be  specified,  as  well  as  the  wall  geometry  and 
boundary  conditions.  The  flow-field  variables  and  boundary  values  are 
updated  point -by -point  from  iteration  to  iteration  until  convergence  is 
reached. 


6.0  RESULTS  AND  DISCUSSION 


In  this  section,  the  techniques  which  have  been  developed  are 
applied  to  the  computation  of  turbulent  internal  flows.  A series  of 
computations  is  discussed  in  the  order  of  increasing  complexity. 

First,  solutions  for  planar  diffuser  flows,  with  constant  eddy  viscos- 
ity, are  described.  The  solution  of  a conical  diffuser  flow,  with 
algebraic  eddy  viscosity  models  and  a wall  matching  procedure,  is 
then  presented.  Because  of  limitations  in  the  algebraic  eddy  viscosity 
models,  solutions  obtained  with  the  two-equation  turbulence  model  and 
the  wall  matching  procedure  are  then  presented.  Finally,  the  two- 
equation  turbulence  model,  along  with  the  numerical  technique  for 
computing  the  wall  region,  is  applied  to  planar  channel  and  planar 
diffuser  flows. 


6.1  SOLUTION  FOR  A PLANAR  DIFFUSER  FLOW  WITH  A 

CONSTANT  EDDY  VISCOSITY 

The  formulation  presented  in  the  previous  sections  was  applied 
to  the  calculation  of  flow  in  the  planar  diffuser  shown  in  Fig.  10. 

The  complete  finite  difference  formulation  with  a constant  eddy  vis- 
cosity is  given  in  Appendix  C.  The  exit  to  inlet  area  ratio  is  two, 
and  the  diffuser  half -angle  is  26.  5 deg.  The  nonuniform  inlet  veloc- 
ity is  represented  by  a parabolic  profile  for  a fully  developed  channel 
flow.  The  inlet  station  is  placed  one  inlet  channel  width  ahead  of  the 
diverging  section  so  that  the  upstream  influence  of  the  diverging  sec- 
tion on  the  inlet  profiles  is  negligible.  A long  exit  section  (8  inlet 
channel  heights)  is  used  in  the  calculation  to  provide  enough  distance 
for  the  flow  to  reach  a parallel  flow  condition,  i.  e. , 8 /8  x = 0.  The 
exit  section,  although  not  always  existing  in  diffuser  applications,  is 
convenient  for  computational  purposes  because  a parallel  flow  bound- 
ary condition  can  be  easily  incorporated  in  the  program. 
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Six  cases  were  calculated  which  correspond  to  the  Reynolds  num- 
bers, based  on  the  inlet  centerline  velocity  (Uc),  the  channel  width 
(h),  and  { v + vt)  of  6,  60,  120,  600,  and  3,000.  The  variation  of  the 
Reynolds  number  was  easily  achieved  by  changing  the  value  of  the 
eddy  viscosity.  The  solutions  are  also  applicable  to  laminar  flows  if 
( v + v^)  is  taken  to  be  v . 

The  calculated  centerline  velocity  distribution  (normalized  by  the 
the  average  inlet  velocity  (U)  is  shown  in  Fig.  11.  At  low  Reynolds 
number,  the  centerline  velocity  remains  at  a constant  value  of  1.  5 in 
the  inlet  section  except  in  the  region  near  the  diverging  section.  The 
velocity  begins  to  decrease  in  this  region  before  the  flow  enters  the 
diverging  section.  The  upstream  effect  is  related  to  the  elliptic 
nature  of  the  flow.  The  disturbance  created  by  the  turning  of  the 
flow  in  the  diverging  section  propagates  upstream  to  cause  the  flow  to 
slow  down.  The  magnitude  of  the  "upstream  effect"  decreases  as  the 
Reynolds  number  increases.  At  a Reynolds  number  of  3,000,  the 
centerline  velocity  varies  only  slightly  from  the  inlet  to  the  exit.  The 
physical  mechanism  of  the  changing  behavior  of  the  flow  in  the  diverg- 
ing section  as  a function  of  Reynolds  number  can  be  understood  more 
easily  by  examining  the  velocity  profiles  and  streamline  patterns 
given  in  Figs.  12.  At  low  Reynolds  number,  the  streamline  follows 
the  wall  smoothly  (Fig.  12a).  At  Re  = 60,  the  flow  separates  from 
the  wall  at  x/h  * 1. 5.  In  the  separated  flow  region,  the  reverse  flow 
velocity  is  very  low  (Fig.  12b).  The  point  of  separation  moves  up- 
stream toward  the  inlet  corner  of  the  diffuser  as  the  Reynolds  number 
increases  (Fig.  12c).  It  is  interesting  to  note  that,  once  the  separation 
point  moves  to  the  upstream  corner,  the  central  flow  becomes  jet-like. 

The  centerline  and  the  wall  static  pressure  distributions  are 
shown  in  Figs.  13  and  14,  respectively.  Direct  comparison  of  the 
two  sets  of  data  shows  very  little  vertical  pressure  gradient.  The 
linear  pressure  drop  in  the  constant  area  inlet  section  corresponds 
to  fully  developed  channel  flow  (Ref.  35).  Since  the  flow  near  the  wall, 
in  general,  turns  more  than  the  flow  near  the  centerline,  the  pressure 
distribution  across  the  diffuser  shows  a larger  nonuniformity  at  the 
inlet  corner  section  than  at  any  other  section.  At  low  Reynolds  num- 
ber, the  pressure  rises  sharply  in  the  diverging  section,  but  the 
rise  tends  to  diminish  as  the  Reynolds  number  increases.  For  the 
jet  flow  cases.  Re  = 3,000  for  example,  the  pressure  remains  con- 
stant which  is  typical  of  completely  "stalled"  diffuser  flows. 

To  investigate  the  quality  and  accuracy  of  the  present  numerical 
formulation,  a conventional  upwind  difference  formulation  was  also 
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employed  in  the  solution  of  the  planar  diffuser  problem.  The  six  test 
cases  were  calculated  by  using  the  upwind  difference  scheme  to  400 
iterations.  The  final  residue  ratio  is  well  below  1 x lO-?.  The  six 
cases  were  re-calculated  by  using  the  upwind  scheme  for  the  200 
iterations,  followed  by  the  present  formulations  to  400  iterations. 

A typical  centerline  velocity  distribution  is  shown  in  Fig.  15.  The 
upwind  difference  scheme  predicts  a faster  velocity  decay  than  the 
present  method.  Typical  behavior  of  the  two  solution  techniques  is 
indicated  in  Fig.  16  by  the  history  of  the  centerline  velocity  conver- 
gence at  x/h  = 2.  0.  The  difference  between  the  two  solutions  becomes 
most  pronounced  at  Re  = 120.  As  would  be  expected,  the  difference 
gradually  diminishes  with  increasing  Reynolds  number.  At  low 
Reynolds  numbers,  i.  e.  , Re  -*■  0,  the  difference  between  the  two 
schemes  also  approaches  zero  because  the  contribution  from  the  con- 
vection terms  vanishes  and  the  flow  is  completely  controlled  by  the 
diffusion  mechanism. 

The  convergence  should  be  considered  in  the  evaluation  of  any 
iteration  process.  The  velocity  residue  ratio  (AU/U)  is  shown  in 
Fig.  17  in  terms  of  the  iteration  number.  The  upwind  scheme  was 
used  for  the  first  200  iterations.  The  residue  ratio  reached  a value 
below  1 x 10"6  in  150  iterations.  The  change  in  numerical  method  to 
the  present  scheme  with  decay  functions  at  the  200th  iteration  causes 
the  residue  ratio  to  rise.  But  it  decreases  rapidly  within  the  next  100 
iterations.  The  solution  obtained  at  the  end  of  400  iterations  is  essen- 
tially the  converged  solution.  Once  the  Q-ijj  solution  is  obtained,  the 
pressure  equation  (Eq.  (13))  is  solved  iteratively.  The  rate  of  con- 
vergence of  the  centerline  pressure  at  x/h  = 2 is  shown  in  Fig.  18. 

One  possible  explanation  for  the  slow  rate  of  convergence  compared 
to  that  shown  in  Fig.  17  is  that  the  boundary  condition  for  the  pressure 
equation  is  of  the  gradient  type  which  usually  requires  many  iterations 
for  convergence. 


6.2  NUMERICAL  SOLUTION  OF  A CONICAL  DIFFUSER 

The  calculations  made  with  a constant  viscosity  model  have  been 
useful  to  define  the  nature  of  the  flow-field  solutions  and  to  gain 
experience  with  the  numerical  method.  However,  real  turbulent 
diffuser  flows  are  poorly  described  by  a constant  eddy  viscosity;  in 
a real. flow,  the  eddy  viscosity  vanishes  as  the  wall  is  approached. 

In  this  section,  algebraic  models  are  used  in  the  calculation  of  the 
separated  flow  in  a 23 -deg  conical  diffuser.  Two  eddy  viscosity 
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models  are  used,  namely  the  Prandtl's  mixing  length  model  and  a 
"convective"  model.  The  region  very  close  to  the  diffuser  wall  is 
analyzed  by  the  use  of  the  law  of  the  wall;  therefore,  the  boundary 
condition  for  the  numerical  solution  of  the  remainder  of  the  flow  field 
is  an  effective  slip  velocity  distribution  along  the  wall.  The  slip 
velocity  is  related  to  the  local  skin  friction  coefficient  by  the  law  of 
the  wall,  so  one  can  assume  either  a skin  friction  coefficient  distribu- 
tion or  a slip  velocity  distribution. 

The  geometry  of  the  diffuser  and  the  computational  grid  are 
shown  in  Fig.  19.  A nonuniform  grid  system  is  used  to  provide 
adequate  flow  definition  near  the  wall.  The  numerical  formulation 
for  a nonuniform  grid  system  does  not  pose  any  serious  difficulties, 
except  that  simple  approximations  for  the  decay  functions,  such  as 
given  by  Eq.  (84),  are  not  readily  available.  Therefore,  the  decay 
functions  must  be  evaluated  at  every  grid  point  with  exponential 
functions.  The  derivation  of  the  decay  functions  for  a nonuniform 
grid  system  is  given  in  Appendix  B. 

The  mixing  length  model  used  in  the  present  analysis  is 


where  S.  = 0. 4 y in  the  region  near  the  wall  and  SL  = 0.09  6 for  the 
remainder  of  the  flow.  The  thickness  (6)  is  an  assumed  function  of 
x.  The  wall  slip  velocity  was  specified  so  that  separation  was  fixed 
at  x/D  = 0.  619.  Predicted  wall  and  centerline  pressure  distributions 
are  shown  in  Fig.  20.  The  wall  pressure  drop  in  the  diffuser  inlet 
section  is  realistic,  but  the  pressure  distribution  downstream  of 
separation  is  unrealistic  in  that  the  wall  pressure  reaches  a maxi- 
mum and  then  decreases.  This  unrealistic  behavior  of  the  pressure 
distribution  is  attributed  to  the  eddy  viscosity  model.  The  mixing 
length  model,  with  a rather  conventional  specification  of  the  length 
scale  distribution,  cannot  provide  an  adequate  solution  for  the  flow 
downstream  of  separation. 

More  reasonable  flow-field  solutions  can  be  obtained  by  specify- 
ing a "convective"  model,  which  approximately  includes  the  tendency 
of  the  viscosity  to  be  convected  along  streamlines  (Fig.  21).  Wall  and 
centerline  pressure  distributions  computed  with  the  convective  eddy 
viscosity  model  are  shown  in  Fig.  22;  the  unrealistic  peak  in  the 
pressure  distributions  has  been  eliminated.  The  predicted  distribu- 
tion of  the  axial  component  of  velocity  is  compared  with  experiment 
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(Ref.  36)  in  Fig.  23.  A comparison  of  the  predicted  and  experimental 
wall  pressure  distributions  is  shown  in  Fig.  24,  along  with  the  assumed 
distribution  of  the  wall  slip  velocity.  The  experimental  results  are 
reasonably  well  predicted. 

In  the  results  which  have  been  discussed,  the  wall  slip  velocity 
distribution  was  held  constant.  The  feasibility  of  relaxing  the  slip 
velocity  distribution  as  the  flow-field  solution  is  relaxed  was  investi- 
gated. A smooth  converged  wall  slip  velocity  distribution  can  indeed 
be  obtained  as  long  as  the  separation  point  is  specified.  And,  in 
principle,  one  can  iterate  on  the  location  of  the  separation  point.  But 
in  practice,  the  logarithmic  wall  region  profile  becomes  singular  at 
the  separation  point.  Perhaps  another  wall  region  profile,  such  as  a 
power  law,  could  be  used  in  the  vicinity  of  separation.  But  the  validity 
of  any  assumed  profile  in  the  vicinity  of  separation  is  questionable. 
Clearly,  an  adequate  solution  requires  that  the  equations  of  motion  for 
the  wall  sublayer  be  solved  along  with  the  solution  of  the  remainder  of 
the  flow  field.  A complete  numerical  solution  requires  a low  Reynolds 
number  turbulence  model  and  coordinate  stretching  in  the  sublayer. 

The  results  presented  in  this  section,  although  not  predictions  in 
the  true  sense  of  the  word,  show  that  accurate  flow-field  solutions 
can  be  obtained  if  reasonable  eddy  viscosity  and  wall  slip  velocity 
distributions  are  assumed.  But  adequate  predictions  of  the  entire 
flow  field  requires  a turbulence  model  which  takes  into  consideration 
that  the  turbulence  is  convected  by  the  mean  flow  field. 

6.3  NUMERICAL  SOLUTION  WITH  A HIGH  REYNOLDS  NUMBER 

TWO-EQUATION  k-e  MODEL  AND  A WALL  MATCHING  PROCEDURE 

In  this  section,  the  two-equation  k-e  model  is  used  for  the  flow 
outside  the  wall  sublayer,  while  the  wall  slip  velocity  boundary  con- 
dition is  retained.  The  purpose  of  this  study  was  to  validate  the 
k-e  model  for  the  high  Reynolds  number  portion  of  the  flow. 

6.3.1  Numerical  Solution  for  a Fully  Developed  Channel  Flow 
with  a Wall  Matching  Procedure 

The  detailed  numerical  formulation  for  the  channel  flow  is  given 
in  Appendix  D.  It  is  commonly  known  that  a region  described  by  the 
law  of  the  wall  exists  in  a turbulent  boundary  layer  when  the  Reynolds 
number  is  high.  The  velocity  profile  is  described  by  the  simple 
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logarithmic  law  (Eq.  (26)).  The  existence  of  the  "law  of  the  wall 
region"  provides  a birdge  between  the  sublayer  and  the  numerically 
calculated  core  region.  The  basic  criterion  which  allows  the  use  of 
Eq.  (26)  in  the  commonly  valid  region  is  that  both  the  velocity  (u) 
and  the  friction  velocity  (v*)  must  be  continuous  at  the  matching 
location  (yQ).  Since  Eq.  (26)  contains  two  unknowns  (u  and  v*),  one 
of  them,  v*  in  this  case,  must  be  determined  iteratively  so  that 
continuous  solution  can  be  obtained.  The  determination  of  v#  re- 
quires the  additional  assumption  that  the  shear  stress  is  constant 
from  the  wall  to  the  matching  location.  The  validity  of  this  assump- 
tion and  its  effect  on  the  accuracy  will  be  discussed  below. 

There  are  two  methods  for  determining  the  friction  velocity; 
the  direct  and  the  indirect: 


In  the  direct  shear  stress  determination,  the  friction  velocity 
v*  is  calculated  from 

v*aJwF=-F7T|  (90) 

at  the  matching  location  (yQ).  The  velocity  gradient  (3u/8  y)  and  the 
eddy  viscosity  are  determined  numerically  from  the  core  region 
solution.  This  simple  and  straightforward  method,  unfortunately, 
does  not  provide  accurate  determination  of  the  friction  velocity  be- 
cause both  the  eddy  viscosity  (v  t)  and  the  velocity  gradient  must  be 
determined  iteratively.  As  a result,  the  value  of  the  friction  veloc- 
ity (v*)  fluctuates  from  iteration  to  iteration.  Hence,  the  velocity 
determined  from  Eq.  (26)  does  not  give  a smooth  convergent  value. 
The  other  factor  which  affects  the  accuracy  of  the  numerically  deter- 
mined velocity  gradient  (8u/8y)  is  the  grid  size.  Coarse  grid  sizes 
can  have  a significant  effect  on  the  numerical  value  of  8u/8  y be- 
cause large  velocity  gradients  normally  exist  in  the  matching  region. 
Thus,  it  is  better  to  replace  the  direct  method  with  the  indirect 
method. 


In  the  indirect  friction  velocity  calculation,  the  method  uses  two 
grid  points  in  the  matching  region  to  determine  the  correct  match- 
ing condition  for  u and  v*.  First,  Eq.  (26)  is  applied  to  a point  at 
(yQ  + 5y)  to  determine  v*  iteratively  as 


(91) 
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where  n is  the  iteration  number.  With  v*  determined,  the  correct  veloc- 
ity at  the  matching  location  is  calculated  from 

U0  s tf*{2.5W9-tf7i>)  + 5.5j  Jo:  specified.  (92) 

a*  do  * ° 

This  method  does  not  require  the  evaluation  of  the  eddy  viscosity  (vt) 
nor  the  velocity  gradient  and  thus  eliminates  the  unnecessary  fluctua- 
tion in  determining  v*  and  uQ.  The  method  has  been  successfully 
employed  to  obtain  convergent  solutions  for  the  channel  flow  and  the 
8-deg  conical  diffuser  flow  described  in  Section  6.3. 2. 

In  addition  to  the  boundary  conditions  for  uQ  and  v *,  one  needs 
to  have  boundary  values  for  the  turbulent  kinetic  energy  and  its 
dissipation.  To  determine  the  turbulent  kinetic  energy  (k)  at  the 
matching  location,  it  is  usually  assumed  that  the  production  of  the 
turbulent  kinetic  energy  is  balanced  with  the  dissipation,  i.  e.  , 

*•(£>*■  e «» 

By  using  Eq.  (93)  and  the  Prandtl-Kolmogorov  eddy  viscosity  formula 
(Eq.  (23)),  one  obtains 

Similarly,  by  combining  Eq.  (93),  the  law  of  the  wall  (Eq.  (26)),  and 
the  definition  of  the  shear  stress  (r),  one  can  determine  eQ  as 


t - j).  fill.')  = ) = ir*  v* 

ce  1 u } z — 


The  numerical  solution  for  a fully  developed  channel  flow  was 
obtained  with  a standard  Gauss -Seidel  iteration  procedure.  At  the 
end  of  each  iteration,  the  boundary  condition  at  the  matching  loca- 
tion is  updated  by  using  Eqs.  (91),  (92),  (94),  and  (95).  The  initial 
guess  is  provided  by  the  logarithmic  law  of  the  wall.  The  initial 
viscosity  is  calculated  from  the  Prandtl's  mixing  length  theory 
(Eq.  (16)).  The  distribution  of  the  mixing  length  is  obtained  from 


l 
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(96) 
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which  is  the  Nikuradse's  formula  for  pipe  flow.  The  initial  turbulent 
kinetic  energy  and  its  dissipation  are  calculated  from 


k = 4 •<#>//£ 


(97) 


and 

c =(!  Jsl 

t > vt 

The  computation  was  carried  out  with  51  grid  points. 


(98) 


In  the  iteration  procedure,  the  eddy  viscosity  was  held  constant 
for  the  first  700  iterations  so  that  the  distributions  of  the  turbulent 
kinetic  energy  and  its  dissipation  can  be  brought  to  a convergent 
state.  The  Prandtl-Kolmogorov  eddy  viscosity  model  is  used  there- 
after with  an  under-relaxation  factor  (rj ) of  0. 1,  i.  e. , 


(99) 


The  convergence  of  the  vorticity  at  the  matching  location  is  shown 
in  Fig.  25  for  a Reynolds  number  of  30,  533  based  on  the  centerline 
velocity  (Uc)  and  the  half  channel  width  (h).  The  solution  proceeded 
through  1, 200  iterations.  The  fairly  constant  value  of  the  vorticity 
indicates  that  the  solution  has  reached  a steady  value.  The  calculated 
velocity  profile  is  shown  in  Fig.  26.  Agreement  with  Laufer's  ex- 
perimental data  (Ref.  37)  at  Re  = 30,  800  is  excellent.  The  match- 
ing location  specified  in  the  calculation  was  yQ  = 0.  06.  The  calculated 
velocity  profile  near  the  wall  follows  closely  the  law  of  the  wall.  A 
small  deviation  from  the  law  of  the  wall  profile  occurs  near  the  cen- 
terline so  that  the  symmetry  condition  can  be  satisfied.  The  log- 
arithmic law  of  the  wall  does  not  satisfy  the  symmetry  condition. 

The  calculated  turbulent  kinetic  energy  distribution  across  the 
channel  is  shown  in  Fig.  27.  The  maximum  TKE  occurs  at  the  match- 
ing location  where  the  profile  closely  follows  the  experimental  data  by 
Clark  (Ref.  38).  The  present  solution  slightly  overpredicted  the  tur- 
bulent kinetic  energy  near  the  centerline.  In  general,  however,  the 
agreement  is  good. 

The  calculated  shear  stress  distribution  is  given  in  Fig.  28, 

The  profile  is  fairly  linear  as  it  should  be  for  a fully  developed 
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channel  flow.  For  a linear  shear  stress  distribution,  the  shear 
stress  at  the  matching  location  is  always  lower  than  that  at  the  wall. 
The  difference  can  be  minimized  by  making  the  location  closer  to  the 
wall.  In  general,  this  can  be  achieved  only  by  increasing  the  Reynolds 
number  of  the  flow  since  the  matching  location  must  be  outside  the 
sublayer  in  order  to  use  the  law  of  the  wall. 

6.3.2  Numerical  Solution  of  an  8-deg  Conical  Diffuser  Flow 
with  a Fully  Developed  Inlet  Velocity  Profile 

One  of  the  few  experimental  investigations  of  diffuser  flow  field 
is  reported  by  Okwoubi  and  Azad  (Ref.  39).  Velocity  distributions 
were  obtained  at  various  stations  in  an  8-deg  conical  diffuser.  The 
inlet  condition  of  the  experiment  is  a well-defined  fully  developed 
pipe  flow  profile,  which  can  be  easily  represented  by  the  present 
numerical  procedure.  In  the  present  numerical  calculation,  the 
inlet  profile  is  specified  through  the  use  of  the  Van  Driest  formula. 

The  inlet  condition  was  allowed  to  relax  to  the  self-consistent  fully 
developed  pipe  flow  profile  in  the  iteration  process.  A coordinate 
transformation  is  necessary  to  map  the  diverging  section  and  the  tail 
section  into  a constant  diameter  pipe.  The  calculation  is  then  per- 
formed in  the  transformed  plane.  The  detailed  numerical  formulation 
of  the  governing  equations  is  given  in  Appendix  E.  The  diffuser 
geometry  is  shown  in  Fig.  29. 

The  calculated  velocity  field  in  the  8-deg  conical  diffuser  is 
shown  in  Fig.  30.  The  development  of  the  velocity  profile  from  that 
of  a fully  developed  channel  flow  into  a free  shear  profile  is  clearly 
demonstrated  in  the  figure.  The  change  in  the  turbulent  flow 
behavior  can  also  be  seen  in  Fig.  31  in  terms  of  the  turbulent 
kinetic  energy  distributions  at  the  inlet  and  a far  downstream  station. 
At  the  inlet  station,  the  TKE  has  a maximum  value  at  the  core- 
sublayer matching  location  near  the  wall  and  monotonically  decreases 
to  the  centerline..  On  the  other  hand,  the  maximum  TKE  appears  in 
the  middle  of  the  diffuser  at  a far  downstream  station  which  is  typical 
of  results  associated  with  a free  shear  profile.  In  Fig.  32,  the 
turbulent  shear  stress  distribution  is  also  given  at  the  two  stations. 
The  linear  shear  stress  distribution  at  the  inlet  gradually  changes 
into  a sine  function  profile  of  a free  shear  layer  at  a far  downstream 
station.  The  location  of  the  maximum  shear  stress  in  the  diffuser 
diverging  section  is  shown  in  Fig.  33.  The  predicted  position  of 
maximum  shear  stress  agrees  fairly  well  with  the  experimental  data 
of  Ref.  39  indicating  that  the  turbulence  transport  is  well  modeled  in 
the  present  numerical  solution.  In  Fig.  34,  the  centerline  velocity 
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distribution  in  the  diverging  section  of  the  diffuser  is  also  shown  to 
agree  well  with  experimental  data.  The  predicted  value,  in  general, 
is  slightly  lower  than  the  experimental  value,  but  the  experimental 
data  show  some  scattering  especially  in  the  middle  of  the  diffuser. 

A comparison  of  the  experimental  and  predicted  velocity  profile 
across  the  diffuser  is  shown  at  four  stations  in  Fig.  35.  The  agree- 
ment between  the  experiment  and  theory  at  the  inlet  and  at  the 
x/D  = 5.95  is  excellent.  At  the  intermediate  stations,  however,  the 
predicted  centerline  velocity  is  about  3 percent  lower  than  the  data. 
This  discrepancy  was  probably  caused  by  the  wall  matching  procedure 
near  the  wall.  The  experimental  data  show  that  the  velocity  profile 
near  the  wall  in  the  diverging  sections  is  not  represented  by  the  log- 
arithmic law  of  the  wall.  In  the  present  approach,  since  the  stream 
function  is  used,  a slight  overprediction  of  the  wall  slip  velocity  in 
the  diverging  section  can  be  magnified  through  the  (1/r)  factor 
associated  with  the  axisymmetric  flow  formulation  (Appendix  E). 

Calculations  for  the  8 -deg  conical  diffuser  flow  were  performed 
with  41  lateral  grid  points  across  the  flow,  which  is  considered  more 
than  adequate  to  provide  accurate  numerical  solution.  The  program 
developed  also  can  be  applied  to  variable  shape  2-D  or  axisymmetric 
channel  flow  problems. 


6.4  NUMERICAL  SOLUTION  WITH  A LOW  REYNOLDS  NUMBER 
TWO-EQUATION  k-e  MODEL 


The  k-e  model,  along  with  the  wall  matching  procedure,  has  been 
shown  to  yield  reasonable  results  for  flows  in  which  the  law  of  the  wall 
is  applicable.  However,  as  pointed  out  in  Section  6.2,  the  validity  of 
the  law  of  the  wall  is  questionable  in  the  vicinity  of  the  separation 
point.  In  order  to  avoid  the  problems  associated  with  the  law  of  the 
wall,  the  whole  flow  field  including  the  sublayer  region  is  solved  by 
the  finite  difference  formulation  so  that  the  point  of  separation  and 
the  separated  flow  field  can  be  predicted.  The  numerical  formulation 
of  the  whole  flow  field  requires  the  use  of  the  low  Reynolds  number 
version  of  the  k-e  model  as  well  as  a sublayer  coordinate  stretching 
technique. 

6.4.1  Numerical  Results  for  a Fully  Developed  Channel  Flow 

The  numerical  formulation  of  a fully  developed  channel  flow  with 
a low  Reynolds  number  two -equation  k-e  model  requires  the  sublayer 
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coordinate  stretching  to  provide  adequate  resolution  through  the  flow 
field.  The  detailed  finite  difference  formulation  of  a fully  developed 
channel  flow  is  given  in  Appendix  F. 

The  velocity  profiles  with  Reynolds  numbers  (Uh/i/)  ranging  from 
1,700  to  207,000  are  shown  in  Fig.  36.  In  the  sublayer  region,  all 
velocity  profiles  follow  the  linear  velocity  distribution,  i.  e. , U + = y+. 
As  the  Reynolds  number  increases,  the  velocity  profile  gradually 
approaches  the  law  of  the  wall  in  the  core  region.  The  velocity  pro- 
files in  the  physical  coordinate  are  presented  in  Fig.  37.  The  profile 
shape  changes  from  near  parabolic  at  the  low  Reynolds  number  to  a 
fuller  profile  at  higher  Reynolds  number.  At  Reynolds  number  1, 657, 
the  profile  agrees  well  with  experimental  data  of  Patel  and  Head 
(Ref.  40).  Results  at  higher  Reynolds  number  also  agree  well  with 
Laufer's  data  (Ref.  37). 

The  turbulent  shear  stress  distribution  is  given  in  Fig.  38.  At 
high  Reynolds  number,  the  turbulent  shear  stress  distribution 
approaches  the  linear  profile  of  the  total  shear  stress  because  the 
contribution  from  the  molecular  viscosity  is  small.  The  maximum 
turbulent  shear  stress  also  appears  near  the  wall.  As  the  Reynolds 
number  decreases,  the  location  of  the  maximum  turbulent  shear 
moves  away  from  the  wall  and  deviates  from  the  linear  distribution 
of  the  total  shear  stress  because  of  the  significant  contribution  from 
the  molecular  viscosity.  The  present  calculated  shear  stress  distribu- 
tion at  Re  = 5,360  agrees  well  with  Eckelmann's  experimental  data 
(Ref.  41)  at  Re  = 5,600, 

The  TKE  distribution  is  given  in  Fig.  39.  The  location  of  the 
maximum  TKE  also  moves  toward  the  wall  as  the  Reynolds  number 
increases.  The  numerical  result  closely  follows  the  experimental 
data  of  Clark  (Ref.  38)  except  near  the  centerline. 

The  effect  of  the  total  number  of  grid  points  on  the  accuracy  of 
the  numerical  solution  is  shown  in  Fig.  40.  With  only  41  grid  points, 
the  calculated  total  shear  stress  deviates  substantially  from  the 
exact  solution.  On  the  other  hand,  the  101 -point  case  is  in  excellent 
agreement  with  the  exact  solution. 

In  Fig.  41,  the  ratio  of  the  centerline  velocity  and  the  mean 
velocity  (averaged  across  the  channel)  is  shown  versus  the  Reynolds 
number.  The  value  of  the  velocity  ratio  (Uc/U)  is  directly  related 
to  the  flatness  of  the  velocity  profile.  For  a uniform  velocity  profile, 
Uc/U  is  equal  to  one.  For  Reynolds  number  below  1,000,  Uc/U  is 
equal  to  1.5,  which  indicates  that  the  velocity  profile  is  parabolic. 
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The  velocity  ratio  decreases  from  a value  of  1.  5 as  the  Reynolds  num- 
ber increases.  The  velocity  ratio  will  presumably  reach  1.0  asymptotic- 
ally as  the  Reynolds  number  goes  to  infinity.  The  skin  friction  coefficient 
is  presented  in  Fig.  42.  The  agreement  between  the  numerical  result  and 
the  experimental  data  (Ref.  42)  in  both  Figs.  41  and  42  is  very  good. 

When  the  low  Reynolds  number  model  is  used  to  calculate  the 
eddy  viscosity  in  the  numerical  iteration,  it  is  necessary  to  use  an 
under-relaxation  factor  on  v ^ as  in  Eq.  (99)  to  provide  smooth  con- 
vergence. The  effect  of  the  relaxation  factor  (r?)on  the  solution 
convergence  is  shown  in  Fig.  43.  Large  oscillations  occur  when  the 
relaxation  factor  (77)  is  greater  than  0.05,  It  was  found  that  the  best 
value  of  77  is  the  inverse  of  the  number  of  lateral  grid  points  which 
in  this  case  is  0.01.  When  the  relaxation  factor  is  greater  than  0.01, 
the  sublayer  thickness  oscillates  as  the  iterations  proceed.  The 
oscillation  phenomenon  becomes  more  pronounced  at  high  Reynolds 
numbers.  Therefore,  the  proper  selection  of  the  under -relaxation 
factor  for  the  eddy  viscosity  associated  with  the  low  Reynolds  num- 
ber two- equation  model  is  very  important  to  ensure  the  smooth 
convergence  of  the  numerical  iteration  process. 

6.4.2  Numerical  Solution  of  Separated  and  Non-Separated  Diffuser 
Flow  in  a Stretched  Coordinate  System 

A series  of  preliminary  calculation  was  made  for  a family  of 
planar  diffuser  flows  to  illustrate  the  nature  of  the  solution.  The 
test  case  selected  is  a two-dimensional  planar  diffuser  with  a 4:1 
aspect  ratio  (Fig.  44)  investigated  by  Reneau,  et  al. , (Ref.  1). 

The  Reynolds  number  based  on  the  inlet  mean  velocity  (Uj)  and  the 
height  (hj)  is  1. 2 x 105.  The  inlet  profile  was  a fully  developed 
channel  flow  profile.  It  was  shown  that  the  maximum  pressure 
recovery  occurred  with  a total  diffuser  angle  of  20  deg.  This  experi- 
mental evidence  is  used  to  provide  limited  verification  for  the  present 
numerical  solutions  because  more  detailed  experimental  data  are  not 
available.  The  coordinate  transformations  are  given  in  Eqs.  (53) 
through  (57),  and  the  set  of  governing  equations  is  given  in  Eqs.  (66) 
through  (70).  Fifty-one  lateral  grid  points  are  used  across  the 
diffuser  which  includes  both  the  core  region  and  the  sublayer  region. 

The  number  of  grid  points  is  considered  adequate  to  provide  a quali- 
tative description  of  the  flow  field  but  not  necessarily  an  accurate 
result.  The  numerical  procedure  is  given  in  Fig.  9.  In  the  itera- 
tion process,  the  first  600  iterations  are  used  to  determine  the  flow 
field  at  the  first  three  stations  so  that  a fully  developed  channel  flow 
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profile  can  be  obtained.  With  the  inlet  condition  known,  the  next  600 
iterations  are  used  to  compute  the  diffuser  flow  field.  The  calculated 
skin  friction  coefficient  is  shown  in  Fig.  45  for  six  different  angles. 
Based  on  the  skin  friction  distributions,  flow  separation  does  not  occur 
in  the  first  three  cases,  namely,  2 Q = 3.  58,  7.15,  and  14.25  deg. 

The  point  of  separation,  which  appears  in  the  last  three  cases,  moves 
upstream  as  the  diffuser  angle  increases.  The  rise  in  the  skin  fric- 
tion coefficient  around  the  inlet  corner  (X/hj  = 2.0)  is  a result  of  the 
elliptic  nature  of  the  flow,  i.  e.  , disturbances  propagate  upstream 
from  the  diverging  section.  Such  behavior  cannot  be  obtained  from 
a conventional  boundary-layer  formulation  in  which  there  is  no  feed- 
back mechanism.  The  skin  friction  coefficient  near  the  exit  corner 
(X/hj  = 6)  in  the  separated  region  shows  some  oscillation  (Figs.  45e 
and  f),  which  indicates  that  the  solution  is  not  fully  converged.  How- 
ever, the  solution  is  stable  upstream  of  the  separation  point. 

The  development  of  the  velocity  profile  in  the  diffuser  is  shown 
in  Fig.  46.  As  the  diffuser  angle  increases,  the  fully  developed 
channel  flow  profile  gradually  develops  into  a wake  profile  near  the 
wall.  For  the  separated  profiles,  the  sublayer  is  so  thin  that  it  looks 
as  if  the  velocity  profile  has  a discontinuity  near  the  wall.  The  veloc- 
ity profile  in  the  sublayer  with  an  enlarged  scale  is  shown  in  Fig.  47 
for  2 0 = 34.7  deg.  The  reverse  flow  velocity  near  the  wall  at 
X/hj  = 6 is  about  14  percent  of  the  local  centerline  velocity.  The 
predicted  turbulent  kinetic  energy  distribution  is  shown  in  Fig.  48. 

The  distinctive  feature  of  the  turbulent  kinetic  energy  distribution 
is  that  the  location  of  the  maximum  TKE  moves  away  from  the  wall 
in  the  diverging  section.  In  general,  the  magnitude  of  the  maximum 
TKE  increases  as  the  diffuser  diverging  angle  increases.  For  2 6 
= 34.7  deg,  the  maximum  TKE  at  X/hj  = 6 is  roughly  double  the 
maximum  inlet  value. 

Although  there  are  no  detailed  flow-field  data  available  to 
verify  the  predicted  flow-field  structure,  the  centerline  velocity 
distribution  is  an  indication  of  the  pressure  recovery.  In  Fig.  49, 
the  centerline  velocity  distribution  is  given  in  terms  of  the  total 
diffuser  angle.  The  minimum  centerline  velocity  at  X/hj  = 6 
occurs  at  2 0 = 20  deg,  which  corresponds  to  the  experimentally 
defined  optimum  diffuser  angle  by  Reneau,  et  al.  (Ref.  1).  A 
comparison  between  the  present  fully  developed  inlet  velocity 
profile  and  that  of  curve  5 in  Ref.  1 is  given  in  Fig.  50. 

To  provide  some  indication  of  the  accuracy  of  the  present 
solution,  the  total  shear  stress  distribution  at  the  inlet  station  is 
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shown  in  Fig.  51.  Although  the  linearity  of  the  total  shear  stress  is 
preserved  in  the  core  region,  the  slope  of  the  profile  is  somewhat 
lower  than  the  exact  analytical  solution.  The  accuracy  check  is 
consistent  with  that  given  in  Fig.  40.  Clearly,  the  number  of  grid 
points  used  in  the  computation  is  only  marginal.  An  increase  to  an 
81-  or  101 -point  system  should  provide  excellent  result  as  was 
found  for  fully  developed  channel  flow  (Section  6.4.1  and  Fig.  40). 


7.0  CONCLUDING  REMARKS 


A numerical  method  has  been  developed  to  provide  the  detailed 
flow-field  structure  of  two-dimensional,  turbulent,  imcompressible 
stalled  and  non-stalled  subsonic  diffuser  flows  with  nonuniform  inlet 
conditions.  The  general  formulation  is  also  applicable  to  a wide 
variety  of  incompressible  internal  flow  problems.  An  important 
feature  of  the  numerical  method  is  the  use  of  the  ''decay  function" 
technique,  in  which  local  analytical  information  is  used  in  the  finite 
difference  formulation  to  ensure  stability  of  the  solution.  A co- 
ordinate transformation  including  sublayer  stretching  was  developed 
to  provide  adequate  flow  definition  throughout  the  whole  flow  field. 
Further  study  of  the  coordinate  transformation  could  provide  opti- 
mization of  the  grid  network. 

A hierarchy  of  solutions  were  obtained,  based  on  turbulent 
transport  models  of  increasing  complexity,  i.  e. , constant  viscosity, 
algebraic,  and  two-equation  models.  In  general,  treatment  of  the 
wall  layer  with  a matching  procedure,  based  on  the  law  of  the  wall, 
yields  fairly  good  results  for  non-separated  flows,  but  is  unsatis- 
factory for  separated  flows.  Therefore,  a method  was  developed  to 
compute  the  whole  flow  field  numerically,  including  the  viscous 
sublayer.  The  approach  involves  coordinate  stretching  and  a low 
Reynolds  number  two -equation  turbulent  model.  Predictions  of 
fully  developed  channel  flows,  obtained  with  the  sublayer  stretching, 
are  in  good  agreement  with  experimental  results.  Predictions  of 
the  performance  of  planar  diffusers  are  also  in  reasonable  agree- 
ment with  experimental  results.  However,  additional  correlation 
between  the  numerical  method  and  experiment  is  needed,  particular- 
ly for  stalled  axisymmetric  flows.  Unfortunately,  the  currently 
available  data  are  very  limited,  particularly  for  evaluating  the 
quality  of  the  flow  structure. 
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The  present  numerical  analysis  should  be  extended  to  include 
compressible  and  non-adiabatic  flows. 
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y+  e yv*/v 

Figure  2.  Turbulent  kinetic  energy  distribution  and  models  for  the 
nonisotropic  part  of  turbulent  kinetic  energy  dissipation 
near  a wall. 
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a.  Uniform  system,  equal  spacing 


b.  Uniform  system,  unequal  spacing 


I 


c.  Nonuniform  system 

Figure  3.  Uniform  and  nonuniform  coordinate  systems. 
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Figure  4.  Body-aligned  coordinate  transformation. 
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Final  Transformed  Domain 
Figure  7.  A complete  transformation  for  a diffuser. 
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Grid  Reynolds  Number,  R 
Figure  8.  Decay  function. 
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Figure  9.  Flow  chart  of  the  numerical  solution  procedure. 
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Axial  Location,  X/h 

Figure  10.  Two-dimensional  planar  diffuser  and  computational  grid  system. 
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a.  Re  = 6 

Figure  12.  Velocity  and  stream  function  distribution  in  a planar 
diffuser. 
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Axial  Location,  X/h 


b.  Re  = 60 
Figure  12.  Continued. 
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c.  Re  = 600 
Figure  12.  Concluded. 


72 


Pressure  Coefficient 


Figure  13.  Centerline  pressure  distribution  in  a planar  diffuser. 
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Pressure  Coefficient 


Axial  Location,  X/h 

Figure  14.  Wall  pressure  distribution  in  a planar  diffuser. 
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Vel 


Figure  15.  Centerline  velocity  distribution  in  a planar  diffuser. 
Re  = 120. 
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Figure  16.  History  of  the  centerline  velocity  convergence 
at  X/h  = 2.0. 
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assure  residue  ratio  at  X/h  = 2.0. 
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Figure  19.  23-deg  conical  diffuser  and  computational  ip-id  in  radial 
direction  (nonuniform  system). 
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Static  Pressure,  p/p 


Figure  20.  Wall  and  centerline  pressure  distribution  with  a mixing 
length  model. 
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Figure  21.  Schematic  of  a convective  eddy  viscosity  model. 
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Static  Pressure,  p/pj 


Figure  22.  Wall  and  centerline  pressure  distribution  with  a convective  model. 


AEDC-TR-76-159 


Radius 


0 1.0  0 1.0  0 0 1.0  0 1.0  0 1.0 

u/uc 

Figure  23.  Velocity  distribution  in  a 23-deg  conical  diffuser  with 
a convective  eddy  viscosity  model. 
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Pressure  Coefficient,  C-  Slip  Velocity 
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a.  Assumed  slip  velocity,  0.1  in.  from  wall 


Axial  Location,  in. 


b.  Wall  pressure  distribution 

Figure  24.  Wall  pressure  and  slip  velocity  distributions  in  a 
separated  conical  diffuser. 
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Vortlcity  at  Matching  Location,  Q(y0) 


Iteration  Number,  n 

Figure  25.  The  convergence  of  the  iteration  process  for  a fully 
developed  channel  flow  calculation. 
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Figure  26.  Velocity  profile  in  a fully  developed  channel  flow, 
with  k-e  model  and  wall  matching. 
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Turbulent  Kinetic  Energy,  k+ 
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Figure  27.  Turbulent  kinetic  energy  distribution  in  a fully 
developed  channel  flow. 
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Figure  28.  Shear  stress  distribution  in  a fully  developed  channel 
flow. 
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Figure  30.  Velocity  distribution  in  an  8-deg  conical  diffuser 
at  various  axial  stations. 
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Turbulent  Kinetic  Energy,  k/v 
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Figure  31.  Turbulent  kinetic  energy  distribution  in  an  8-deg 
conical  diffuser. 


91 


Centerline  Velocity, 


Figure  34.  Centerline  velocity  distribution  in  an  8-deg  conical  diffuser. 
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Radial  Location,  r/0 
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a.  X/D,  = -0.595 

Figure  35.  Comparison  between  the  present  calculated  velocity 
distribution  and  experimental  data. 


95 


Radial  Location,  r/0.5  Dj 


AEDC-TR-76-159 


b.  X/D,  = 0.4 
Figure  35.  Continued 
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c.  X/D,  = 4.05 
Figure  35.  Continued. 
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Figure  37.  Velocity  profiles  in  a fully  developed  channel  flow. 
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Figure  38.  Turbulent  shear  stress  distributions  in  a fully 
developed  channel  flow. 
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Figure  40.  Effect  of  the  total  number  of  grid  points  on  the  total 

shear  stress  distribution  in  a fully  developed  channel  flow. 
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Re  = u*  h/v 


Figure  41.  Comparison  between  the  calculated  and  experimental  centerline 
velocities  in  a fully  developed  channel  flow. 


Skin  Friction  Coefficien 
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Figure  42.  Comparison  between  the  calculated  and  experimental  skin 
friction  coefficient  in  a fully  developed  channel  flow. 
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Wall  Vorticity,  ft  x 10" 


Number  of  Iterations,  n 

Figure  43.  Effect  of  viscosity  relaxation  on  the  convergence  of  the 
wall  vorticity. 
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Figure  44.  A 2-D  planar  diffuser  with  a 4:1  aspect  ratio. 
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Skin  Friction  Coefficient 
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a.  2 6 - 3.58  deg 

Figure  45.  Skin  friction  coefficient  in  a 2-D  diffuser  flow. 
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Figure  45.  Continued. 
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c.  2 0 = 14.25  deg 
Figure  45.  Continued. 
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X/hj 

d.  2 e - 21.2  deg 
Figure  45.  Continued. 


e.  2 6 = 28.1  deg 
Figure  45.  Continued. 
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Figure  47.  Sublayer  velocity  distribution,  (2  0=  34.7  deg). 
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Figure  48.  Turbulent  kinetic  energy  distribution  in  a 2-D  diffuser. 
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r/0.5  h (X) 

b.  2 6 = 7.15  deg 
Figure  48.  Continued. 
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c.  2 6 = 14.25  deg 
Figure  48.  Continued. 
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r/0.5  h(X) 

d.  2 0 = 21.2  deg 
Figure  48.  Continued. 
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e.  29  = 28.1  deg 
Figure  48.  Continued. 
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f.  26  = 34.7  deg 
Figure  48.  Concluded. 
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Figure  51.  Accuracy  of  the  inlet  profile  based  on  the  total  shear 
stress  distribution. 
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Table  1.  Comparison  between  Two  Low  Reynolds  Number  Models 
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APPENDIX  A 

DERIVATION  OF  A COORDINATE  TRANSFORMATION 
WITH  A SUBLAYER  STRETCHING 


The  purpose  of  the  transformation  is  twofold:  (1)  to  provide  ade- 
quate resolution  in  the  sublayer  region  which  normally  occupies  only 
about  1 percent  of  the  flow  field,  and  (2)  to  obtain  a smooth  stretching 
of  the  coordinates  from  the  sublayer  to  the  core  region.  Because  of 
the  complexity  of  the  velocity  profile  near  the  wall,  separate  trans- 
formation functions  are  derived  for  the  sublayer  and  the  core  region. 
They  are  matched  smoothly  at  a proper  location  (yQ)  so  that  the 
functions  are  continuous  up  to  the  second  derivatives. 

TRANSFORMATION  FUNCTION  IN  THE  SUBLAYER 


Inside  the  sublayer,  the  velocity  distribution  is  linear.  A tangent 
function  is  a suitable  stretching  function  for  the  sublayer  region: 

y tan  (A-l) 

The  stretching  factor  (9y/9y),  which  must  be  adjusted  to  provide  good 
resolution  of  the  sublayer  velocity  profile,  is  related  to  the  slope  of 
the  velocity  profile  by 


jj  I _ ^ 9HJ  or*2 


(A-2) 


The  condition  determines  the  coefficient  (a)  in  terms  of  v*2/j/.  The 
coefficient  ((3)  will  be  determined  through  the  matching  procedure  at 


y0- 


TRANSFORMATION  FUNCTION  IN  THE  CORE  REGION 

In  the  core  region,  there  is  no  need  to  have  a large  coordinate 
stretching  because  the  profile  is  rather  smooth.  Therefore,  the 
stretching  factor  (dy/dy)  must  gradually  decrease  to  unity  from  the 


131 


AEDC-TR-76-159 


sublayer  to  the  centerline.  In  the  present  analysis,  a hyperbolic 
function  is  selected  for  this  purpose,  i.  e. , 

|i=  c + tank  (A-3) 

Upon  integration,  Eq.  (A-3)  becomes 

* ln(cosh(%-tjt))  4?  (A-4) 

Equation  (A-4)  is  the  transformation  function  for  the  core  region. 


THE  MATCHING  CONDITION 


The  matching  conditions  are  provided  at  the  location  (yQ)  by  the 
continuity  of  the  function,  the  first  and  the  second  derivatives.  With 
the  continuity  of  the  second  derivatives  at  yQ,  i.  e. , 


core 


2f 


sublayer 

% 


(A-5) 


the  coefficient  (/3)  can  be  determined  from  Eqs.  (A-l),  (A-4),  and 
(A-5)  as 


2«>tp  sec*p^ -tan  = o 


(A-6) 


For  the  continuity  of  the  first  derivatives,  the  condition  becomes 


aT  «i|s“Miyer 
’5  h “ *1 


(A -7) 


Equation  (A -7)  determines  the  coefficient  (C)  as 


c = <*  sec4(a& 


(A-8) 
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Finally,  the  continuity  of  the  transformation  functions  provides  the 
following  relation  for  the  constant  (F),  i.  e. , 

F = -C%  +(-2otanpJ.  (A-9) 

By  setting  y = 1.0,  the  transformed  centerline  location  (ymax)  can 
be  determined  from  Eq.  (A-4)  as 

C?m«-  *■  fC0Sh( WS. »*F  “>  <A-10> 

Equations  (A-2),  (A-6),  (A-8),  and  (A-9)  permit  the  transformation 
functions  (A-l)  and  (A-4)  to  be  written  in  terms  of  v^/v  and  yQ, 
which  are  characteristic  parameters  of  the  sublayer. 


NUMERICAL  PROCEDURE 


The  numerical  procedure  used  to  calculate  the  coordinate  trans- 
formation functions  is  outlined  in  Fig.  A-l.  First  the  parameter 
(v*2/i/)  is  determined  based  on  some  flow  condition,  such  as  the 
inlet  condition  to  a diffuser.  The  stretched  sublayer  thickness  is 
then  set  equal  to  unity  so  that  enough  resolution  can  be  obtained.  In 
addition,  the  number  of  grid  points  is  specified  in  the  transformed 
y coordinate.  The  coefficients  (a,  C,  and  F)  can  be  calculated 
directly  from  Eqs.  (A-2),  (A-8),  and  (A-9),  respectively.  On  the 
other  hand,  coefficients  such  as  £S  and  ymax  must  be  determined 
iteratively  by  Newton's  method,  i.  e. , 


Prwi  58  ft. 


f(gsL 

f'w 


(A-ll) 


where  0n+1  represents  the  (n+l)th  value  and  0n  represents  the  old 
value. 


The  functions  (f  and  f')  are  obtained  from  Eq.  (A-6)  as 

2<*pn  Sec2p„  -tan 1 3n  -I  <A-i2) 

^Pb)  sfpl  * 2o<  SeciPn  +PnSecpr)+2flttan*p„) 

^ (A-13) 

where  yG  =1.0  has  been  used  in  the  derivation.  The  converging 
solution  for  (3  is  obtained  by  applying  Eq.  (A-ll)  successively  from 
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an  initial  guess  until  a convergence  criterion  is  reached.  Since  the 
Newton's  method  converges  fairly  rapidly  with  a good  initial  guess, 
the  maximum  number  of  iteration  is  often  specified  to  terminate  the 
calculation. 

The  coefficient  (F)  can  also  be  determined  in  a similar  manner, 
i.  e.  , 

~ ~ l ( ?«««,») 

W’nt'  ‘ 

where  f and  f'are  now  derived  from  Eq.  (A- 10)  as 

^ * “ * + ^ ^max,io  + ^ (COSh  C </max,n  “ 

f ” C + "tunh  ( ^m<>n  “ • ) 

With  all  the  coefficients  determined,  one  can  proceed  to  determine 
the  uniform  grid  spacing  (Ay  = (ymax/JNM)),  the  stretched  coordinate 
(y  = (N.  Ay)),  and  the  physical  coordinate  (y)  by  Eqs,  (A-l)  and  (A-4). 
The  coordinate  transformation  as  a function  of  v*2  jv  is  shown  in 
Fig.  A-2. 


(A-14) 

))+p  (A-15) 

(A-16) 


The  transformation  factors  derived  from  Eqs.  (A-l)  and  (A-4) 
are  in  the  sublayer,  0 < y < y^: 


Q/( {p- -*/•*■■ 


(A-17) 


and  in  the  core  region,  yG  < y < y 


max' 


| 

W'  C + tanh(J-5#) 

%t'  / = - [ Sf 1 ^ 1 f H. ' >]  /[ C * tank  (Zj-i)] 
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Figure  A-1.  Numerical  procedure  to  determine  the 
coordinate  transformation  function. 
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Figure  A-2.  Coordinate  transformation  as  a function  of  the 
sublayer  parameter  (v*2/u). 
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DERIVATION  OF  DECAY  FUNCTIONS 
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The  derivation  of  decay  functions  is  illustrated  by  using  the  steady- 
state  vorticity  equation.  The  derivation  of  the  decay  function  for  the 
time -dependent  equation  is  given  in  Ref.  30.  The  one -dimensional 
vorticity  equation  is 


7**0  <B_1) 

The  finite  difference  expressions  with  the  decay  functions  for  a uniform 
grid  system  are 


^jfi  “ + -Qj-i  | 

w ’% 


iQ-  „ j-n  ~ Ai-i  l 

n ^ zSy  Fj 

The  resulting  finite  difference  equation  is 


(B-2) 


J / V \ nS*\~& 

<*j  1 y O'  2sy 


■j-f 


(B-3) 


The  function  (Fj)  is  not  explicitly  used  in  Eq.  (B-3)  because  only  one 
decay  function  is  required  for  an  equation  with  only  two  terms.  For 
a nonuniform  grid  system,  one  obtains  the  finite  difference  equation. 

f flfH-A-yfr,  - (A-AnVf  »>  i fV)  1 _ 0 

I (sy,  + sy,)/2  % (fa,+ss.i  J 

where  6yj  and  6y2  are  the  grid  spacing  between  (j+1)  and  (j),  (j)  and 
(j-1)  grid  points,  respectively.  After  re-arrangement,  Eq.  (B-4) 
becomes 


(V/.+SHz) 


2 (f  faty t)  J 


(B-5) 


The  local  analytical  solution  of  Eq. 
ing  that  v is  locally  constant,  i.  e. , 
Eq.  (B-l)  is 


(B-l)  can  be  obtained  by  assum- 
v..  The  analytical  solution  of 

J 


a s c,  + ca-  e 


(B-6) 
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The  local  boundary  condition  is 

^ “ 0 » jfl « 


(B-7) 

The  final  expression  for  y = 6y2  is  obtained  by  combining  Eqs.  (B-6) 
and  (B-7). 

, Aw  i 

Qj  = flj.,  + j (b-8) 


The  decay  function  (GJ  can  be  easily  obtained  by  direct  comparison 
of  Eqs.  (B-5)  and  (B-B).  The  result  is 


r c*»,+s8t)  r sa»  AW . 

I * <■>!  MM  1 

J ( 


(B-9) 


for  a uniform  grid  system,  6y^  = 5y2  = 5y,  and  Eq.  (B-9)  reduces  to 

2(e<vis» 


g.'-L—L.IKH  -i)  ) 

J 0f$**  t - 1 ) J 

or  in  terms  of  the  grid  Reynolds  number  (Rj  = (v/v)j  6y), 


(B-10) 


2(e^J-/) 
e2l?J‘  -i 


) 


(B-l  1) 


Equations  (B-9)  and  (B-ll)  represent  the  decay  functions  in  a non- 
uniform  and  a uniform  grid  system,  respectively. 
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The  governing  equations  for  a two-dimensional  planar  steady- 
state  flow  of  an  incompressible  fluid  with  a constant  viscosity  can 
be  written  as: 


flu  av 

fly 


= o 


(C-l) 


1A 


9X  Yfly 


f flX  3Xl  fly*' 


fl*U 


Vu 


v fltj  J ay  r 1 3x*  ay* 1 


(C-2) 

(C-3) 


Equations  (C-l)  through  (C-3)  can  be  written  in  terms  of  the  vorticity- 
stream  function  formulation  as 


£“♦*1  --L(  iiM+vH)  = 0 

ax*  df  y '•"ax  S')' 


*X£, 

A second-order  pressure  equation  can  also  be  derived  as 


ax*  ay* 


2 


(C-4) 
(C  — 5) 


(C-6) 


where  the  source  term  of  Eq.  (C-6)  has  been  written  in  terms  of  the 
derivatives  in  the  y- direction  for  ease  in  the  numerical  computation. 
The  velocity  field  (u,v)  is  obtained  from  the  relation 


U = 


it 


(C-7) 


The  corresponding  finite  difference  equations  for  the  general  network 
shown  in  Fig.  C-l  are: 
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Vortichy  Equation  (ft) 


(fiitlj  + flKj)  1 ,U.i 

si*  (C_8) 

+ Gj  *lj  ***  ' 


where  the  decay  functions  Gj  and  Gj  are  defined  as 

<rL  = (2/rl)(i-  2(e^-l)/(eR- 1 

Qj  = (2/Rj)(j-2  (eRj-i)/fe2Pj-0) 
s ( u/p  )£  j Sx  > K j 3 ( v/v  )i,j s i 


(C-9) 

(C-10) 

(C-ll) 


Stream  Function  Equation  {\p) 

(W2,ry + V**1  * Wi.jH-W.i'Vi.j-'  V*s’  • --i.j 

(C-12) 


i4s-n. 


Pressure  Equation  (p) 


Velocity  Relation  (u,  v) 


(C-14) 


When  Eqs.  <C-8),  (C-12),  and  (C-13)  are  solved  iteratively  by  Gauss- 
Seidel  iterative  method,  the  corresponding  successive  substitution 
formula  can  be  derived  as 
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. (Cr$i*i.i  + +c3*fr,j*i  SOURCE) 

s cu 

(C-15) 

The  corresponding  C^,  C2,  C3,  C4,  SOURCE,  and  CU  terms  for 
each  function  are: 


ft-Equation 

1 

C,  =■  i/et  - Ri  A 

C2  b i/gl  + fJi/i 

C3  = ( l/Gj  - Rj/a  )•(  S*AS  )2 

C4  .(l/e-j  + R//2  )•(«/*»)* 

SOURCE  * 0 

CU  » (2/S£)  + (2/Sj) 

and  the  decay  functions  are  approximated  as 


= |.0-0.062S-(Rt)a 

> IRiKz 

2 1 

, IRiU2 

IRJ  " (R L)‘ 

* 1.0  - 0.0625- fRj)2 

> lRjl<2 

s 2 1 

, JRjUa 

~ IRjl  “(Rj)1 

^-Equation: 


(C-16) 


(C-17) 


c,  =1  c4=(ex/sy)2 

C2  = 1 c*u  = 24  2 (Z*M) 

c3  = LSt/sy  )2  SOURCE  * ilg  CSX)2  (C-18) 
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p-Equation: 


C, 

C2  »l 

C3  » (Sx/sy)1 


CA  » (8x/sy)2 

CU  a 2+2  (Sx/gy)2 


SOURCE  = iffUi.j  4 )/**•]. 


The  boundary  conditions  used  for  the  planar  diffuser  calculation  are: 


UPSTREAM  BOUNDARY  CONDITION 

A fully  developed  parabolic  velocity  profile  is  specified  at  the 
inlet  of  the  diffuser  in  the  present  analysis  to  represent  the  non- 
uniform  inlet  condition.  However,  any  other  inlet  profile  can  also 
be  specified.  The  corresponding  vorticity  and  the  stream  function 
distributions  are  derived  from  the  velocity  profile.  The  inlet  static 
pressure  profile  is  assumed  to  be  uniform. 

DOWNSTREAM  BOUNDARY  CONDITIONS 


When  the  length  of  the  exit  section  of  the  diffuser  is  long  enough 
for  a parallel  flow  to  be  established,  the  parallel  flow  condition  can 
be  expressed  as 


9A  _ W 
9x  " dx  m 


(C-20 


The  corresponding  downstream  condition  for  the  pressure  can  be 
derived  from  the  momentum  Eq.  (A-2)  as 


lts^A 

9X 


(C-21) 
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SYMMETRY  CONDITION 


Along  the  line  of  symmetry,  both  the  vorticity  and  the  stream 
function  are  set  equal  to  zero. 


The  pressure  boundary  condition  along  the  line  of  symmetry  can 
be  obtained  from  Eq.  (C-3)  by  integration  in  y direction,  i.  e. , 


(C-22) 


Subscripts  1 and  2 represent  the  centerline  and  its  neighboring  grid 
point,  respectively  (see. Fig.  C-l).  By  assuming  that  the  integrands 
are  linear  function  of  y,  one  obtains 


(C-23) 


where  the  continuity  equation  has  been  used  for  dv/a y.  The  second 
and  the  third  terms  inside  the  bracket  represent  the  higher  order 
term  which  vanishes  in  the  limits  as  the  ay  goes  to  zero. 


WALL  BOUNDARY  CONDITION 

At  the  wall,  the  velocity  components  vanish  and  the  stream  func- 
tion is  a constant  value  along  the  solid  wall,  i.  e.  , 

u = v » o , 

The  vorticity  at  the  wall  can  be  obtained  by  the  integration 

rW  -W  ... 

where  w and  p represent  the  wall  and  its  neighboring  points, 
that  the  integration  can  be  performed  in  either  the  x-  or  y-direction. 
Assuming  that  the  integrand  in  Eq.  (C-25)  is  linear,  one  obtains 

+ (c-26) 

Similarly,  the  wall  boundary  condition  can  be  derived  by  integration 
of  the  momentum  Eqs.  (C-2)  and  (C-3)  in  either  the  x-  or  y-direc- 
tion. For  example,  the  pressure  boundary  condition  along  a 


(C-24) 

(C-25) 

Note 
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diverging  diffuser  wall  is  obtained  by  integrating  Eq.  (C-2)  in  the 
x-direction. 


VSX 

2 


f/M)  +/KL)  1 + 

I**' ty  7»J + 2 

[UW-V2“1 


gx 

2 


(C-27) 


For  the  horizontal  inlet  and  exit  walls,  the  condition  (Eq.  (C-27)) 
can  be  reduced  to 


K - 


T>SX 


<c-28> 


The  formulation  is  fairly  general  in  that  it  can  be  used  for  calculat- 
ing most  2-D  planar  flows  with  proper  assignment  of  the  boundary 
conditions. 
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APPENDIX  D 

FORMULATION  OF  A FULLY  DEVELOPED  CHANNEL  FLOW  WITH  A HIGH 
REYNOLDS  NUMBER  TWO-EQUATION  k-6  MODEL  AND  WALL  MATCHING 


For  high  Reynolds  number  fully  developed  channel  flow,  the 
governing  equation  is  one  dimensional.  Since  the  sublayer  is  replaced 
by  a law  of  the  wall  matching  procedure,  the  effect  of  the  molecular 
viscosity  can  be  ignored.  The  complete  set  of  governing  equations 
is: 


Vorticity  Equation  (£2) 


rii  JLj 
ar*  H 1 *r 


AslM  + *.iLi±l>LAl 

r J dr  * vt  1 r ar  rl  J 

+ ' /.»!&.*:*  i =, 
Pt  l ar*  ar  J 


(D-l) 


Stream  Function  Equation  (^) 


ar1  Str  ar 


+ r si  * o 


(D-2) 


Velocity  Recovery  (u) 


U 


Turbulent  Kinetic  Energy  Equation  (k) 


1 r > or 


(D-3) 


ar 


>rl  Vtl  r Jar  'ar'  ^ 


(D-4) 


Turbulent  Kinetic  Energy  Dissipation  Equation  (e) 


»r*  Ptl  ar'n  r‘«j  J»r 


(D-5) 


P ran dti- Kolmogorov  Eddy  Viscosity  Model 


* - 9*  T 


(D-6) 
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The  constants  used  are  = 0.09,  Cj  = 1.44,  and  C2  = 1.92.  The 
boundary  conditions  along  the  centerline  can  be  written  as 


¥*■=  constant 


it  sli 
«r  ar 


0 


(D-7) 


The  matching  condition  near  the  wall  is  provided  by  the  law  of  the  wall. 
The  boundary  condition  and  the  related  variables  are 


v*=  Up/(2.5  in  + 

U **  U**(2.5  In  (V*$e/}))  + 5.5) 

“0 


k * v*2/V^T 

£ = v*i/(oAt  yc) 
dr 

(D-8) 


where  yQ  is  the  matching  location  measured  from  the  wall  and  yp  is 
the  distance  one  point  next  to  the  matching  point.  The  governing 
equations  written  in  the  standard  form  are 


+ d a 0 


(D-9) 


where  coefficients  a and  d are  given  in  the  following  table  for  a fully 
developed  channel  flow  high  Reynolds  number  model: 


4> 

a 

d 

n 

1 l 2%  1 f 9*4  3ttl 

r ar  ”rar  arJ 

V 

r8ii 

B 

mu 

e 

B 
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The  general  finite  difference  formulation  of  Eq.  (D-9)  becomes 


l sr2  JGj  j|  2sr  J J 

where  the  decay  function  (Gj)  is  calculated  from 


= 0 


(D-10) 


1.0-  0.0625  (Rj)1 

> iRjl  < 2 

2 1 

|Rj  1 £ 2 

> 

IRjl  ‘ (Rj)2 

(D-ll ) 

aj  • sr 

The  source  terms  (dj)  are  calculated  from  the  conventional  central 
difference  scheme.  Note  that  the  coefficient  (aj  and  Rj)  do  not 
necessarily  represent  the  physical  convective  velocity  and  the 
Reynolds  number.  The  coordinate  system  used  and  the  variation  of 
the  eddy  viscosity  both  contribute  to  the  "effective"  Reynolds 
number  (Rj). 
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APPENDIX  E 

FORMULATION  OF  2-D  OR  AXISYMMETRIC  TURBULENT  FLOWS  WITH  A 
TWO-EQUATION  HIGH  REYNOLDS  NUMBER  TURBULENCE  MODEL 


The  governing  equations  for  2-D  or  axisymmetrical  flows  in 
physical  coordinates  are: 

Vorticity  Equation  (Q) 

1 »r  ' *'  r> 
t J.  I,  ££  . ♦ 20.)  + 2 ( M. . 2“  )l  = 0 

^['3Xa  9r*;i3r  3x9T  1 3T  9K  ' I 


Stream  Function  Equation  (0) 


1 9xl  3rl  J r ®r  " rsi 


(E-l) 


(E-2) 


Velocity  Recovery  (u,  v) 


7®  ar 


t 

v~  ft  ax 


Prandt I- Kolmogorov  Eddy  Viscosity  Model 


H-*t 


Turbulent  Kinetic  Energy  Equation  (k) 


9llt  l x/rn 


(E-3) 


<E-4) 


(E-5) 
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Turbulent  Kinetic  Energy  Dissipation  Equation  (e) 


+ 

where 


{5?+iFrTlC  Cv 

* /-r 


(E-6) 


5 u = 0.09  , Ct  = /.44  , Cx  = /. 92  , =/./ 


For  diffuser  flows,  a transformation  is  applied  to  the  above  set 
of  equations,  i.  e. , 


rmiko 


X ax 


(E-7) 


where  (r,  x)  represents  the  physical  coordinate  system  and  (¥,30 
represents  the  transformed  coordinate  system.  The  chain  rule  and 
the  corresponding  transformation  factors  are 

JL  = J_+ 

9X  3X  3X  l9X  1 


j!_  - i.  , f * i z(™)  + r— i— 

9xl  " a**  1 *x  ' af*  + 2 lax  '3*97  'a? 

_a 

ar  “ 9r  ; 

JL-  » (*L)\2L 
»ra  l«r;af» 

a*  ,ar  n 9*  /2La/2Li  . J_  fUL) 
9xar  B 'ar  Jp”1,  a?  laxar' 


(E-8) 


where  the  transformation  factors  are 


(2L)S-L- 
1 ar ' $(x) 


flL)  = -?  /.^A) 

rijg<xi' 


(E-9) 


axa'~  ' ax  'I  g<xr  “ r *■  $fx>  "(  $rx>  0 


laxar;  ljg(x)»J 
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The  governing  equations  written  in  the  transformed  coordinates  in  the 
standard  form  of: 


are 


(E-10) 


O -Equation 


a,  = i 


tj  = -^  { [ v - * f - f >4  K|f) * [u -2  ( ‘ (|Ji 


° Z(»x'33?3?  * r 1 v+  »?  'ar'  r I %V*xl  'w'lP 

Z(ax>  jx a'p  +t3x*,(3)£  } (ar'  afri  35*  +(axJ?7 

+ 2 / iJL  (2l , + ( 21  )(2t)  £1 1 il  (£L  ))((&  « . « 

Sx3r  9r^  3rM*x  *rl  *r  3x»r VV^r'*?  95? 
-<$>£>}  (E-U) 

-Equation 


b|  = 0 

d =-£-u. +i(|f)^^+rc-A  (E-i2> 
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k-Equation 


Q,  * I 


k ■ i{rv-(#)|^-f  4ic£wu-$*<g>$*& 
-»<£>} 


a,  * i 


<ufu  .[$(*).$.&)£]'} 


-c,% 


<VJc 


* %• 


Velocity  Recovery  (u,  v) 


U = 
* - 
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APPENDIX  F 

FORMULATION  OF  A FULLY  DEVELOPED  CHANNEL  FLOW  WITH  A LOW 
REYNOLDS  NUMBER  TWO-EQUATION  k-e  MODEL  AND  SUBLAYER  COORDINATE 

STRETCHING 


The  complete  governing  equations  for  a fully  developed  channel 
flow  are: 

Vorticity  Equation  (ft) 


ar1 


i £* 
(P+i^)  3rl 


•Jl 


+ $ 


D.  M iyi 

r ar 


0 


(F-l) 


Stream  Function  Equation  (<//) 


Velocity  Recovery  (u) 


*rl 


1 r >r 


u 


.LIE 

j-5  ar 


Prandtl- Kolmogorov  Eddy  Viscosity  Model 


(F-2) 


(F-3) 


(F-4) 


Turbulent  Kinetic  Energy  Equation  (k) 

fk  —i 

*rl“ 


Turbulent  Kinetic  Energy  Dissipation  Equation  (e) 


1 »r  r('  't’jar +(0t4)(ar' 

-iiit  } = 0 

q*  J 


■>r»  l »r 1 rl  si'Jjr 


k *r 


L 

k (P  + i^/e^) 


= 0 


(F-5) 


(F-6) 
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The  coefficients  (C^,  Cj,  C2,  and  ae  ) are 

5u=  A/C3Ol00  + /i/0.2  7)] 

C,  = /.44. 

C2  = f.92  C l.o-  0.3 exp (-^5]  (F"7) 

®e  = 1. 1 

a = <nnr 3 /V  , r h k */()>  e ) 

For  a given  coordinate  transformation,  such  as 

7=7(1-)  (f-8) 


the  above  set  of  equations  can  be  written  in  the  standard  form, 

sty  „ 


1 CL  — 

r ®r 


+ d = 0 


(F-9) 


The  coefficients  (a  and  d)  are  listed  in  the  following  table  for  a fully 
developed  channel  flow  low  Reynolds  number  model: 


4> 

a 

d 

a 

1 # sfc  S (*+>*>] 

(^)l  »f  t F.  J 

-Fz 

XI  / ^ P a^tx 
(v+>£)1  a?1  ha  ®r  J 
. sil  f 1 s$  1 1 

0>  + »tHr97  F,  CF,)1/ 

V 

-c  *s  1 
Fl'f  F. 

rE  » 
CF.)1 

1 

1 [ s?t  S 

*7  »■  Fi  J 

-Fa 

* (w)* 

0>+^)  l*r  ' 

l<  an  1 

91  J(F.)* 

e 

1 ( 

C 9**^  /9U»* 

(P+^/Sf)  1 Jr'V 
S ()>♦*/«)  r 

V+^t/^'37  # 

- C * / » 1 

r F,  Fi 

(P+Jfc/ep1  »t  ;(F0l 

F,  - s?/ar  , Fa  s (^r/arO/Car/ar)1 
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A 

a 


al»  a2 


C 

ju 

C 

D 

d 

e 

F 

Gi.  G, 
h 

JNM 

k 

L 

SL 

n 

P 


NOMENCLATURE 

Parameter 

Coefficient 

Coefficients 

Constant 

Coefficient 

Coefficients 

Skin  friction  coefficient 

Pressure  coefficient 

Eddy  viscosity  coefficient 

Transform  coefficient 

Local  diameter 

Source  term 

Difference  in  finite  difference  equations 
Transform  parameters 
Decay  functions 

Channel  or  2-D  diffuser  height,  y 

Number  of  grid  points  along  the  r or  y coordinate 

Turbulent  kinetic  energy 

Differential  operator 

Prandtl  mixing  length 

Iteration  number 

Pressure 
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R 

Parameter,  k^/ve 

Re 

Reynolds  number 

Ri*  Rj 

Grid  Reynolds  number 

r,  r 

Radial  coordinate  in  physical  and  transformed 
plane,  respectively 

rmax 

Diffuser  wall  coordinate  in  transform  plane 

s(x) 

Diffuser  wall  coordinate  in  physical  plane 

u 

Axial  velocity 

u 

Average  axial  velocity  at  diffuser  inlet 

U,  V,  w 

Velocity  components 

U',  V',  w' 

Trubulent  velocity  components 

u+ 

Sublayer  velocity,  u/v  * 

X 

Axial  coordinate 

a,  3 

Transform  parameters 

7 

Trucation  error 

6 

Index,  zero  for  planar  configurations,  1 for  axisym 
metric  configurations  or  incremental 

p «i> 
0 *'1'* 

Boundary  layer  displacement  thickness 

e 

Isotropic  part  of  the  total  turbulent  energy 
dissipation,  k3^2/i 

6 T 

Total  turbulent  kinetic  energy  dissipation 

f] 

Relaxation  factor 

20 

Total  diffuser  divergence  angle 

V 

Molecular  viscosity 
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Vt 

P 

CTe 

T 

* 

SUBSCRIPTS 


c 

f 

I 

max 

o 

P 

ref 

sep 

T 

w 


Eddy  viscosity 

Density 

Constant 

Shear  stress 

Dependent  variable 

Stream  function 

Vorticity 


Centerline 
Finite  difference 
Inlet  section 
Maximum 

Sublayer,  core  region  solution  matching  location 

Neighboring  point 

Reference 

Separation 

Total 

Wall 
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